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Chapter  1 

Introduction 


Presented  is  a  review  of  the  classical  lattice-gas  method  that  deals  with  an  artificial  many-body  sys¬ 
tem  of  particles  that  has  severely  discretized  microscopic  dynamics  and  that  behaves  like  a  viscous 
Navier-Stokes  fluid  in  the  long  wavelength  hydrodynamic  limit.  We  explain  and  analytically  quan¬ 
tify  how  the  artificial  lattice-gas  system  behaves  and  we  derive  a  set  of  criterion  that  specifies  under 
what  conditions  it  can  serve  as  an  appropriate  model  of  a  viscous  and  compressible  fluid.  Then,  we 
show  how  the  lattice-gas  algorithm  works  using  two  test  models.  Finally,  we  compare  the  numerical 
predictions  obtained  from  a  variety  of  different  simulations  of  these  two  test  models  to  the  respec¬ 
tive  analytical  predictions  we  previously  obtained  for  these  models.  The  resulting  numerical  and 
analytical  predictions  are  in  good  agreement  in  all  cases,  but  this  is  only  after  many  failed  attempts 
that  were  incrementally  corrected  over  time  by  removing  flaws  from  the  derivation  of  the  analytical 
predictions  as  well  as  removing  numerical  bugs  in  the  implementation  of  the  algorithm  and  data  col¬ 
lection  methodology.  Therefore,  the  reason  for  the  consistently  good  agreement  between  numerical 
and  analytical  predictions  is  that  the  derived  criteria  set  has  been  so  sharply  delineated  that  we  now 
know  with  great  accuracy  how  to  initialize  the  numerical  model  within  a  narrow  parameter  regime 
where  the  lattice-gas  system  is  operative.  If  the  initial  state  of  the  lattice-gas  system  is  outside  this 
narrow  operating  regime,  the  numerical  predictions  are  not  at  all  in  agreement  with  the  analytical 
predictions  and  the  behavior  of  the  long  wavelength  modes  in  the  system  can  no  longer  be  classified 
as  hydrodynamical.  We  have  not  attempted  to  catalog  any  of  the  non-hydrodynamical  behaviors  of 
a  classical  lattice-gas  system.  Instead,  we  have  chosen  to  pursue  a  narrower  goal,  which  as  it  turns 
out  is  computationally  more  difficult  to  pursue,  where  we  run  the  algorithm  only  in  a  parameter 
regime  where  it  behaves  like  a  fluid. 


1.1  Background 

Much  information  about  the  details  of  the  microscopic  state  of  a  many-body  system  is  not  relevant 
to  the  hydrodynamic  behavior  of  the  many-body  system  at  the  macroscopic  scale.  So  in  numerically 
simulating  the  macroscopic  scale  fluid  behavior,  it  is  possible  to  arrange  for  the  computer  to  keep 
track  of  many  fewer  particles  than  are  in  the  actual  system.  In  a  typical  lattice-gas  simulation  of  fluid, 
there  are  about  109  particles.  In  the  simplest  type  of  simulation,  all  the  particles  are  indistinguishable 
and  move  independently  of  one  another,  except  that  groups  of  particles  may  collide  together  when 
they  arrive  at  the  same  point.  Over  a  decade  ago,  a  classical  lattice  gas  in  two  spatial  dimensions  was 
found  to  behave  like  a  viscous  Navier-Stokes  fluid  at  the  macroscopic  scale  by  Wolfram  [  ]  and  by 
Frisch,  Hasslacher,  and  Pomeau  [2].  This  is  known  as  the  FHP  lattice  gas.  Soon  after  this  discovery, 
a  classical  lattice  gas  was  found  to  model  three-dimensional  fluids  [3].  In  the  simulation,  there  may 
also  be  fixed  obstacles  with  which  the  particles  have  perfectly  elastic  collisions.  For  example,  one  can 
simulate  vortex  shedding  in  a  fluid  flowing  around  a  fixed  object  [4,  5,  6].  The  value  of  shear  viscosity 
for  classical  lattice-gas  fluids  has  been  studied  by  many  researchers,  including  Wolfram  [1],  Frisch 
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et  al.  [3],  Somers  and  Rem  [7],  Henon  [8].  In  the  case  where  there  is  an  attractive  force  between 
particles  spatially  separated,  multiphase  fluid-like  behavior  is  observed  [9,  10,  11,  12,  13,  14,  15,  16] 
(this  is  a  principal  application  of  the  lattice-gas  method) . 

Microscopic  classical  lattice-gas  models  of  Navier-Stokes  fluids  are  now  well  understood.  The 
lattice-gas  method  has  been  undergoing  improvement,  beginning  in  the  mid  1970’s  up  to  the  present 
day  by  many  researchers,  including  Hardy,  de  Pazzis  and  Pomeau  [17],  Kadanoff,  McNamara,  and 
Zanetti  [18,  19,  20,  21],  Boghosian  et  al.  [22,  23,  24,  25],  Boon  et  al.  [26,  27],  Ernst,  Das,  Brito,  et 
al.  [28,  29,  30,  31],  Henon  [8,  32],  Doolen,  S.  Chen,  et  al.  [33,  34,  35,  36,  37,  38],  Frisch,  Pomeau, 
d’Humieres  et  al.  [2,  3],  Appert,  Zeleski,  and  Rothman  et  al.  [9,  10,  11,  12,  13],  and  Yepez  [14, 
15,  16,  39,  40].  This  is  by  no  means  either  an  exhaustive  list  of  all  researchers  or  publications  in 
this  subject.  An  exhaustive  preprint  archive  is  maintained  at  Los  Alamos  National  Laboratory  (see 
“http:  /  /  xyz.lanl.gov/archive  /  comp-gas” ) . 


1.2  Algorithmic  Scheme 

A  lattice  gas  is  a  system  of  identical  particles  where  the  particles  move  on  a  discrete  spatial 
lattice.  The  spatial  lattice  is  an  array  of  points  (which  is  also  referred  to  as  sites  or  nodes),  arranged 
in  a  regular  crystallographic  fashion,  and  the  lattice  appears  exactly  the  same  from  whichever  of  the 
points  the  array  is  viewed.  This  kind  of  lattice  is  called  a  Bravais  lattice.  In  the  case  of  a  single-speed 
lattice  gas,  as  the  name  implies,  all  the  particles  in  the  system  move  at  the  same  speed:  c  =  -,  where 
I  is  the  distance  between  two  neighboring  sites  of  the  lattice  and  r  is  the  time  it  takes  for  a  particle 
to  hop  from  one  site  to  another  site  in  the  nearest  neighboring  vicinity.  In  a  lattice-gas  system,  all 
particles  move  at  same  time  and  then  they  collide.  The  propagation  step  is  called  the  streaming 
phase.  The  propagation  time,  r,  for  particle  streaming  is  also  the  update  time  for  the  entire  system 
because  we  imagine  that  the  collision  phase  of  the  update  procedure  happens  instantaneously  and 
homogeneously  across  the  entire  system  of  particles.  The  smallest  possible  mean  free  path  length 
for  a  particle  in  a  lattice-gas  fluid  is  on  the  order  of  the  grain  size  of  the  lattice,  l.  There  exists  a 
new  and  unique  global  arrangement  of  particles,  referred  to  as  a  state  of  the  system,  following  each 
and  every  update  time  period  r.1 

The  local  state  of  each  particle  is  specified  by  a  position  coordinate  and  a  momentum  vector.  We 
can  think  of  each  local  state  of  the  system  as  having  a  unique  position  in  the  lattice  and  a  unique 
orientation  or  displacement  vector  (corresponding  to  a  momentum  vector  for  particle’s  motion).  As 
a  particle  moves  through  the  spatial  lattice,  it  hops  from  local  state  to  local  state.  Each  local  state 
is  like  a  container  that  can  temporarily  hold  a  particle  that  is  moving  in  a  particular  direction. 

In  the  simplest  type  of  lattice-gas  model,  no  more  than  one  particle  can  occupy  a  single  local 
state.  However,  more  than  one  particle  can  reside  at  a  single  site  at  any  one  time  since  there 
are  B  number  of  local  states  per  site.  Each  local  state  holds  a  particle  at  that  site  moving  in  a 
unique  direction.  Hence,  the  maximum  number  of  particles  that  can  coexist  on  a  single  site  equals 
the  number  of  nearest  neighboring  sites,  since  it  is  possible  that  particles  can  hop  from  all  nearest 
neighboring  sites  to  one  particular  site  all  at  once.  Of  course,  the  minimum  number  of  particles  that 
can  reside  at  a  site  is  zero,  which  can  happen  when  all  the  particles  at  that  site  move  away  and  no 
new  ones  enter  from  any  of  the  neighboring  sites.  So  the  information  needed  to  encode  a  particle’s 
occupancy  of  a  local  state  is  a  single  classical  bit  associated  with  that  site  and  that  local  state.  If 
the  bit  is  on,  a  particle  is  there.  If  the  bit  is  off,  then  no  particle  occupies  that  local  state. 


1.3  History  of  Lattice-Gas  Developments 

Let  me  briefly  review  some  of  the  developments  of  the  lattice-gas  method  applied  to  fluid  dy¬ 
namics  simulation.  An  overview  of  the  lattice  gas  fluids  has  been  given  by  Boghosian  [41].  Lattice 

1  The  fact  that  every  state  of  the  system  is  unique  follows  from  the  principle  of  detailed-balance  which  is  obeyed 
during  the  on-site  collisions  that  may  occur  after  multiple  particles  arrive  at  a  single  site. 
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gases  are  a  special  case  of  cellular  automata,  originally  introduced  by  von  Neumann  and  Ulam  in 
1948  [  12]  and  popularized  in  the  1980’s  by  Fredkin  [43]  and  by  Stephen  Wolfram  [44,  45].  A  broad 
treatment  of  the  cellular  automata  subject  is  presented  by  Toffoli  and  Margolus  in  their  book  on 
cellular  automata  machines  [46].  Following  the  cellular  automata  paradigm,  lattice  gases  are  suited 
to  fine-grained  parallel  processing,  also  called  massively  parallel  processing. 

A  simple  lattice-gas  model  of  discrete  molecular  dynamics  on  a  square  lattice  was  analytically 
investigated  in  the  early  1970’s  by  the  French,  J.  Hardy,  O.  de  Pazzis,  and  Yves  Pomeau  [17]. 
Computer  implementations  were  not  carried  out  until  a  decade  later.  Their  model,  known  as  the 
HPP  model,  was  the  first  classical  lattice  gas  to  reproduce  hydrodynamic  behavior  at  the  macroscopic 
scale.  In  the  late  1970’s,  cellular  automata  research  was  underway  at  the  Information  Mechanics 
Group  at  MIT  on  reversible  computation  by  Fredkin,  Toffoli,  and  Margolus  [  17,  43,  48].  They  built 
fine-grained  special-purpose  machines  to  simulate  physics-like  models  [49,  48].  A  review  of  the  kind 
of  cellular  automata  modeling  done  in  the  early  1980 ’s  is  given  by  Vichniac  [50].  Stephen  Wolfram,  a 
visitor  at  the  Information  Mechanics  Group  in  1983,  popularized  cellular  automata  as  simple  models 
of  self-organization  amenable  to  statistical  mechanics  analysis  [44,  51  ]. 

In  1985  Wolfram  completed  the  first  hydrodynamic  lattice-gas  simulations  on  a  triangular  lattice 
[52]  on  the  first  Connection  Machine — at  that  time,  lattice  gases  were  a  very  appropriate  application 
for  the  bit  oriented  single  instruction,  multiple  data  Connection  Machine  [53] .  After  visiting  the  MIT 
Information  Mechanics  Group  and  seeing  a  type  of  HPP  lattice-gas  running  on  the  cellular  automata 
machine  CAM-5  of  Toffoli  and  Margolus  [54,  49] ,  Pomeau  was  inspired  by  seeing  hydrodynamic-like 
behavior  (for  example,  the  superposition  of  sound  pulses)  occurring  at  the  macroscopic  scale  by 
simple  rules  homogeneously  applied  at  the  microscopic  scale.  By  1986  Frisch,  Hasslacher,  and 
Pomeau  had  reported  the  existence  of  an  isotropic  two-dimensional  lattice  gas  on  the  triangular 
lattice  [2]  that  is  described  by  the  Navier-Stokes  equation  (1.4).  Their  model  is  referred  to  as 
the  FHP  model.  Accompanying  the  seminal  1986  FHP  paper  was  a  paper  by  Margolus,  Toffoli, 
and  Vichniac  on  cellular  automata  machines  for  fluid  dynamics  modeling  [55].  The  contribution 
of  Margolus  et  al.  was  meant  to  complement  the  work  of  Frisch  et  al.,  pointing  out  the  benefit  of 
dedicated  computational  hardware  for  lattice-gas  models.  In  the  same  year  Wolfram  completed  a 
detailed  treatment  of  the  basic  theory  of  discrete  lattice-gas  fluids  using  symmetry  considerations 
and  worked  out  the  mesoscopic  scale  description  in  the  Boltzmann  approximation  [1]. 

By  1987  the  lattice-gas  methodology  was  extended  to  model  three-dimensional  flows  by  Frisch, 
d’Humieres,  Hasslacher,  Lallemand,  Pomeau  and  Rivet  [3].  The  minimal  lattice  found  was  the 
face  centered  hypercubic  (FCHC)  lattice.  The  FCHC  lattice  is  four-dimensional  with  24-nearest 
neighbors.  It  is  projected  onto  three  dimensions  in  a  simple  fashion  by  limiting  the  depth  of  the 
fourth  dimension  of  the  simulation  volume  to  one  lattice  link.  Much  effort  was  spent  on  finding 
optimal  collisions  to  minimize  the  viscosity  of  the  fluid  [8]. 2 

It  was  realized  in  the  late  1980’s  by  Rivet  and  Frisch  [56]  and  by  McNamara  and  Zanetii  [57] 
that  a  lattice  gas  could  be  simulated  directly  at  the  mesoscopic  scale  using  the  lattice  Boltzmann 
equation.  In  place  of  the  discrete  microscopic  representation,  one  begins  at  the  mesoscopic  scale  by 
representing  the  probability  of  a  particle  occupying  a  local  state.  It  has  the  advantage  of  eliminating 
inherent  noisy  fluctuations  in  the  simulation,  but  at  the  expense  of  discarding  the  particle-particle 
correlations. 

In  1954,  D.  Bhatnager,  E.  Gross,  and  M.  Krook  showed  how  the  collision  integral  in  the  Boltz¬ 
mann  equation  can  be  reduced  to  a  simple  diagonal  form  [58];  this  is  termed  the  BGK  approximation. 
A  cutoff  to  the  collision  integral  is  made  by  neglecting  high  order  particle-particle  correlations  to 
cast  the  Boltzmann  equation  as  an  approximate  partial  differential  equation.  Using  the  BGK  ap¬ 
proximation  of  the  collision  intregal,  Chen,  S.  Chen  and  Mattaeus  showed  how  one  could  remove 
the  anomalies  (such  as  non-Galilean  invariance)  that  occur  at  the  macroscopic  scale  by  tailoring  the 

2  This  task  is  difficult  because  the  FCHC  lattice  gas  has  224  or  16.7  million  local  configurations.  In  practice,  all 
possible  collisions  are  not  included  in  a  simulation  because  of  the  large  demand  for  local  memory  needed  to  pre-store 
all  the  necessary  collisional  events  in  table  look-up  format.  To  ease  memory  loads,  Somers  and  Rem  used  lattice 
isometries  to  reduce  the  size  of  look-up  tables  [7].  An  implementation  of  FCHC  on  the  CAM-8  was  carried  out  in 
1995  by  Alder,  Boghosian,  Flekkoy,  Margolus,  and  Rothman  [22]. 


6 


CHAPTER  1.  INTRODUCTION 


functional  form  of  the  equilibrium  occupancy  probability  [59].  The  drawback  of  the  BGK  approxi¬ 
mation  is  that  it  violates  detailed  balance,  but  it  has  computational  advantages  over  modeling  lattice 
gases  at  the  microscopic  scale  because  of  a  significant  viscosity  reduction.  Subsequently,  a  lattice 
Boltzmann  BGK  model  for  compressible  fluid  dynamics  was  proposed  by  Alexander,  H.  Chen,  S. 
Chen,  and  Doolen  [60]  and  another  lattice  Boltzmann  BGK  model  for  thermohydrodynamics  was 
proposed  by  Alexander,  Chen,  and  Sterling  [61]. 

The  lattice  Boltzmann  equation  in  the  BGK  approximation  has  been  a  useful  contribution  of 
the  lattice  gas  community  to  high-performance  computational  physics  [62]. 3  The  lattice  Boltzmann 
equation  (implemented  within  one  hundred  lines  of  code  in  a  parallel  language  such  as  FORTRAN 
90)  allows  one  to  efficiently  model  fluids  with  near  second-order  convergence  using  an  explicit  time- 
step  scheme.  Martinez,  Matthaeus,  S.  Chen,  and  Montgomery  compared  the  lattice  Boltzmann 
BGK  method  to  the  well  known  spectral  method  and  found  comparable  numerical  efficiency  in  both 
methods  [63].  The  lattice  Boltzmann  method  was  extended  to  the  simulation  of  multiphase  and 
multispecies  applications  by  Shan  [64,  65,  66,  67],  Grunau  [68],  and  Yeomans  et  al.  [69,  70,  71]. 
Because  of  the  practicality  of  the  lattice  Boltzmann  method,  we  have  explored  ways  to  restore 
detailed  balance,  in  particular  by  using  an  unconditionally  stable  collision  process  that  uses  unitary 
matrices  to  produce  the  outgoing  configuration  of  occupancy  probabilities. 

Other  areas  of  lattice  gas  research  include:  thermohydrodynamics  [38,  28,  27,  61,  63],  immiscible 
fluids  [9,  37,  72,  73],  reaction-diffusion  systems  [74,  75,  76],  magnetohydrodynamics  [77,  78,  79, 
80],  flow  through  porous  media  [81,  36],  and  renormalized  kinetic  theory  [82,  83,  84,  85,  86,  87, 
88].  Numerical  measurements  taken  from  classical  lattice  gas  simulations  are  generally  in  excellent 
agreement  with  mean-field  theory  predictions  and,  in  the  rare  instance  when  this  is  not  the  case,  with 
more  exact  field  theoretic  calculations  [82,  85,  88].  Lattice  gases  simulate  physical  systems  while 
keeping  multiparticle  correlations.  The  phenomenon  of  long-time  tails  in  the  velocity  autocorrelation 
function  [89,  90,  91]  has  recently  been  observed  in  lattice  gases  by  KirkPatrick  [82]  and  Brito  [83,  85]. 

A  good  review  of  the  lattice  gas  subject,  with  particular  emphasis  on  interfaces,  phase  transitions, 
and  multiphase  flow,  has  recently  been  presented  by  Rothman  and  Zaleski  [13].  Additionally,  a 
comprehensive  bibliography  of  the  subject  has  been  compiled  by  Doolean  [92]. 


1.4  Lattice-Gas  Applications  to  Multiphase  Fluid  Dynamics 

Microscopic  lattice  gas  models  of  multiphase  fluids  were  introduced  by  Chen  et  al.  [3  ]  and  by 
Appert  and  Zaleski  [10,  11,  12],  and  Yepez  [14,  15,  16].  It  is  known  that  interparticle  potentials 
can  be  modeled  by  including  a  single  anisotropic  nonlocal  interaction  in  the  lattice  gas  dynamics  for 
discrete  momentum  exchange  between  particles.  The  simplest  model  of  this  kind  is  the  Kadanoff- 
Swift-Ising  model  [93] .  A  nonlocal  interaction  of  the  type  was  used  in  a  lattice  gas  scheme  by  Appert 
and  Zaleski  [10]  in  1990  to  cause  an  attractive  force  between  particles  giving  rise  to  an  athermal 
liquid-gas  phase  transition.4 

To  simulate  the  correct  macroscopic  dynamics,  the  interaction  range  for  the  momentum  exchange 
must  be  much  smaller  than  any  scale  related  to  the  interface  region  that  exists  between  the  two  phases 
in  order  for  the  multiphase  lattice  gas  to  be  operative.  For  example,  the  interaction  range  must  be 
much  smaller  than  the  radius  of  curvature  of  a  drop  or  bubble  and  much  smaller  than  the  wavelength 
of  a  capillary  wave  or  gravity  wave  travelling  along  the  interface. 

In  the  bulk  region,  the  lattice-gas  collision  operator  for  the  nonlocal  interactions  vanishes,  and 
only  the  local  collisions  determine  the  hydrodynamic  behavior.  Furthermore,  the  rheology  of  mul¬ 
tiphase  dynamics  is  driven  by  low  Reynolds  number  flows.  The  rheology  of  droplets  (for  example 

3  The  first  lattice  Boltzmann  simulations  for  three  dimensional  flows  with  high  Reynolds  numbers  (about  50,000) 
were  presented  in  June  1993  at  the  International  Conference  on  Pattern  Formation  and  Lattice-Gas  Automata  spon¬ 
sored  by  the  Fields  Institute,  Waterloo,  Canada. 

4  A  nonthermal  lattice-gas  is  one  with  intensive  quantities  for  the  pressure  and  density,  but  no  intensive  quantity 
related  to  temperature.  This  is  because,  a  nonthermal  lattice-gas  is  one  where  all  particles  move  at  a  single  speed  and 
a  particle’s  mass  and  momentum  are  uniquely  defined,  but  its  energy  is  not. 
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diffusive  growth,  coalescence,  or  droplet  breakup  by  high  shear  flows)  is  a  complicated  hydrody¬ 
namic  process  yet  well  approximated  as  a  nonturbulent  one.  Complex  fluid  behavior  of  this  sort 
is  characterized  by  a  well-defined  local  equilibrium  that  exists  throughout  the  entire  fluid  volume. 
The  nonlinear  behaviors  driven  by  convection  or  by  surface  energetics  occurs  while  the  multiphase 
fluid  is  in  local  equilibrium  everywhere  (particularly  in  the  bulk  regions)  and  slowly  approaches  a 
global  thermodynamic  equilibrium.  The  local  equilibrium  of  the  fluid  is  characterized  by  a  nearly 
isotropic  occupancy  of  local  momentum  states.  At  zeroth  order,  the  fluid  is  locally  at  rest  with  the 
probability  of  occupancy  being  uniformly  distributed  in  momentum  space.  The  nonequilibrium  fluid 
is  characterized  by  a  small  first-order  correction  away  from  the  local  equilibrium.  Momentum  slowly 
diffuses  through  the  fluid  system  (quantified  by  the  transport  coefficient  for  viscosity)  in  a  random 
walk  fashion. 

Classical  lattice  gases  use  only  a  single  bit  to  represent  the  occupancy  of  a  particle  in  a  local 
state  whereas  in  molecular  dynamics  codes  a  few  hundred  bits  are  used  (six  floating-point  numbers 
for  position  and  momentum).  To  do  physical  modeling  at  the  microscopic  and  mesoscopic  scales, 
molecular  dynamics  is  the  appropriate  modeling  tool.  However,  at  large  hydrodynamic  scales,  a 
classical  lattice  gas  is  the  modeling  tool  of  choice.  Predicting  the  hydrodynamics  behaviors  of 
microemulsions  (three  species  fluid  of  water,  oil,  and  a  surfactant  with  dipolar  molecules)  have 
demonstrated  that  lattice  gases  are  more  efficient  that  molecular  dynamics.  No  competing  high 
level  scheme  is  known  for  microemulsions.  The  microemulsion  lattice  gases  have  been  studied  by 
Boghosian,  Coveney,  and  Emerton  [24].  They  have  mapped  out  the  transition  between  phases  as  a 
function  of  surfactant  concentration  and  domain  growth  laws,  which  have  been  found  to  be  in  good 
agreement  with  experiment  [25].  Lattice  Boltzmann  simulations  of  microemulsions  has  been  done 
by  Yeomans  et  al.  [94] . 


1.5  Limitations  and  Drawbacks  of  Classical  Lattice  Gases 

Far  away  from  equilibrium  where  local  equilibrium  conditions  on  the  particle  configurations 
are  violated,  classical  lattice  gases  behave  in  unphysical  ways.  These  behaviors  are  not  signs  of 
instabilities,  but  indicate  that  far  away  from  equilibrium  artifacts  caused  by  the  discreteness  of  the 
microscopic  dynamics  can  arise  at  the  macroscopic  scale. 

The  classical  lattice  gas  is  like  its  classical  molecular  dynamics  counterpart.  The  available  number 
of  particles,  per  computer  simulation,  is  still  far  too  few  (on  the  order  of  billions)  compared  with 
the  vast  numbers  of  particles  in  any  natural  situation  (on  the  order  of  Avagadro’s  number)  it  is 
trying  to  represent.  Classical  lattice  gases  fail  to  adequately  capture  turbulence  within  large-scale 
hydrodynamics  motions  because  of  limitations  of  available  memory  resources  in  classical  computers. 
Orszag  and  Yakhot  have  consequently  argued  that  classical  lattice  gases  are  not  as  efficient  as  the 
competing  high  level  CFD  methods  which  make  more  efficient  use  of  available  memory  [95]. 

As  well  as  limited  in  spatial  and  temporal  resolution,  classical  lattice  gases  possess  noisy  fluctua¬ 
tions  [33] .  These  fluctuations  are  a  mechanism  whereby  the  lattice  gas  explores  different  metastable 
states  [13].  They  have  the  negative  aspect  of  effectively  reducing  the  classical  simulation’s  macro¬ 
scopic  scale.  If  the  system  is  ergotic5,  we  have  one  of  two  choices  to  remove  the  noise  in  any 
measurement.  We  can  either  increase  the  spatial  size  of  the  lattice  to  allow  for  more  spatial  and 
temporal  coarse-grain  averaging  or  increase  the  number  of  sample  runs  with  different  initial  con¬ 
ditions.  The  latter  is  a  means  of  ensemble  averaging.  In  either  case,  significant  computational 
resources  must  be  expended  to  remove  noisy  fluctuations  instead  of  expending  these  resources  on 
increasing  the  resolution  of  the  simulation. 

Some  possibilities  have  been  explored  to  try  to  avoid  the  noisy  fluctuations  in  classical  lattice  gas 
simulations.  As  mentioned  above,  in  lattice  Boltzmann  simulations,  noisy  fluctuations  in  the  system 

5  A  lattice  gas  is  said  to  be  ergotic  if  numerical  results  obtained  by  course-grain  averaging  over  a  single  large 
microscopic  realization  are  identical  to  the  numerical  results  obtained  by  ensemble  averaging  over  many  microscopic 
realizations. 
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have  smaller  sizes  than  in  a  system  governed  by  microscopic  lattice-gas  rules.  Yet  the  lattice  Boltz¬ 
mann  method  suffers  from  numerical  instabilities  typically  encountered  in  finite-difference  methods. 
The  reason  for  this  comes  about  by  the  method’s  lack  of  detailed  balance,  or  even  its  lack  of  the 
weaker  condition  of  semi-detailed  balance,  in  its  BGK  collision  operator.6  Since  it  is  essentially  a 
first-order  finite-difference  method  [63],  the  explicit  lattice  Boltzmann  BGK  method  is  not  more 
efficient  than  state-of-the-art  implicit  and  second-order  convergent  computational  fluid  dynamics 
methods,  for  example  methodes  employing  either  multiscaling  or  curvilinear  adaptive  meshing. 


1.6  The  Navier-Stokes  Equation 

The  long  wavelength  hydrodynamic  behavior  of  a  many-body  system  of  particles  can  be  modeled 
at  the  macroscopic  scale  by  an  effective  field  theory,  a  set  of  coupled  partial  differential  equations. 
The  smooth  fields  of  mass  density,  p ,  and  flow  velocity,  v,  obey  a  mass  continuity  equation  and 
a  viscous  Navier-Stokes  fluid  equation  of  motion.  There  is  also  a  parabolic  heat  equation  for  the 
energy  density,  yet  for  simplicity,  I  will  not  consider  the  heat  equation  here,  and  instead  I  shall 
consider  an  athermal  fluid. 

Because  the  mass  increase  within  a  region  1Z  is  entirely  accounted  for  by  the  flux  of  particles  into 
1Z  through  its  boundary  dlZ ,  the  p  and  v  fields  obey  the  continuity  equation 


dtp+di(pvi)=  0.  (1.1) 

This  is  the  first  equation  of  motion.  Here,  the  shorthand  notation  for  partial  derivatives  is  used: 
dt  =  d/dt  and  di  =  d/dxi.  The  field  equation  embodying  Newton’s  second  law,  for  a  region  7 Z 
expressing  the  change  in  the  momentum  density  in  terms  of  the  stress  applied  at  the  boundary  dlZ , 
is  Euler’s  equation 

dt{pvi)  +  djllij  =  0.  (1.2) 

Now  following  Landau  and  Lifshitz  [96],  the  momentum  flux  density  tensor  is  written  as' 

2 

lit-;  =  PSij  +  pviVj  -  r](diVj  +  djVi  -  —dkvk6ij)  -  (Sijdkvk.  (1.3) 

The  viscous  stress  tensor  is  <jb  =  p(diVj  +  djVi  —  jjdkvk5ij)  +  C$ijdkvk,  where  p  and  (  are  the 
transport  coefficients  for  the  shear  viscosity  and  bulk  viscosity,  respectively,  and  D  is  the  number  of 
spatial  dimensions  of  the  system.  The  first  two  terms  in  Equation  (1.3)  represent  the  ideal  part  of 
the  momentum  flux  density  tensor,  which  is  sum  of  the  pressure  term,  P,  plus  the  convective  term, 
pvv,  which  is  nonlinear  in  the  velocity. 

In  general  the  pressure,  P,  is  a  function  of  the  mass  density  field,  p  =  p(x,t),  and  for  a  thermal 
fluid  it  also  is  a  function  of  the  temperature  field,  T  =  T(x,t).  The  pressure  tensor  is  diagonal 
because  the  fluid  is  isotropic.  P  =  P(p ,  T)  is  termed  the  equation  of  state.  For  a  neutral  fluid 
comprised  of  independently  moving  particles,  the  pressure  depends  linearly  on  the  mass  density, 
P  =  c^p,  where  cs  is  the  speed  of  sound  in  the  fluid.  In  a  thermohydrodynamic  system,  the  sound 

speed  is  temperature  dependent,  cs  =  \J (where  kB  is  the  Boltzmann  constant  and  m  is  the 
mass  of  a  single  particle).  In  this  case  the  pressure  obeys  the  well  known  ideal  gas  law ,  P  =  nkBT , 
where  n  =  -&-  is  the  particle  number  density.  For  an  athermal  hydrodynamic  system  (one  where 

6  Mass  and  momentum  are  only  conserved  to  within  the  precision  of  the  floating-point  representation.  If  the  value 
of  the  single-particle  distribution  function  at  some  site  is  close  to  either  zero  or  one,  it  is  possible  that  the  occupancy 
probabilities  there  will  become  either  negative  or  greater  than  one  because  of  numerical  round-off  errors  .  When 
either  of  these  situations  arise,  the  lattice  Boltzmann  simulation  will  eventually  become  unstable  everywhere  and  the 
values  of  the  distribution  function  will  diverge  until  a  numerical  underflow  or  overflow  event  occurs.  Usually  the  BGK 
collision  operator  becomes  unstable  in  a  region  with  a  high  density  gradient,  for  example  at  an  interfacial  boundary, 
or  in  a  region  a  with  a  high  velocity  shear. 

7  For  non-divergent  flow  {djVj  =  0)  in  the  incompressible  fluid  limit,  Equation  (1.3)  is  fl ,,  =  PSij  +  pViVj  + 
p(diVj  +  djVi).  Furthermore,  the  term  rjdtVj  vanishes  in  the  Euler  equation  in  this  limit. 
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the  system  is  at  uniform  homogeneous  temperture,  and  where  heat  transport  is  neglected),  cs  is  a 
constant. 

Substituting  Equation  (1.3)  into  Euler’s  equation  Equation  (1.2),  gives  us  the  second  equation 
of  motion  for  a  viscous  isotropic  fluid 


p  ( dm  +  vjdjVi )  =  -diP  +  pvd2Vi  +  (C  +  ^)  didjvj-  (1-4) 

This  is  the  called  the  Navier-Stokes  equation.  In  Equation  (1.4),  p  is  the  shear  viscosity  and  £  is  the 
bulk  viscosity.  The  transport  coefficient  for  momentum  diffusion,  v  =  ^ ,  is  the  kinematic  viscosity.  It 
gives  a  measure  for  the  rate  of  decay  of  local  shears  in  the  fluid  and  determines  how  fast  a  perturbed 
fluid  will  relax  from  an  anisotropic  flow  profile  at  the  macroscopic  scale  to  an  isotropic  steady  state 
profile.  Both  the  shear  viscosity  and  the  bulk  viscosity  cause  damping  of  compressional  waves  in 
the  mass  density  field.  The  shear  viscosity  alone  causes  damping  of  shear  waves  in  the  momentum 
density  field.  In  general,  for  a  nonisotropic  fluid,  there  may  also  exist  a  cubic  viscosity.  However,  in 
our  case  we  shall  deal  with  isotropic  fluids  where  the  shear  and  cubic  viscosities  coincide. 


1.7  Dimensionless  Numbers 

Let  L  and  T  denote  the  characteristic  length  and  time  scales,  respectively,  of  a  hydrodynamic 
scale  fluctuation.  That  is,  L  and  T  are  quantities  characterizing  the  fluid’s  configuration  at  the 
macroscopic  scale.  Examples  of  the  characteristic  length  scale  for  hydrodynamic  flow  are  the  wave¬ 
length  of  a  compressional  wave  in  the  mass  density  field,  the  wavelength  of  a  shear  wave  in  the 
momentum  density  field,  or  the  diameter  of  a  fluid  vortex.  The  mean  free  path  is  the  average 
distance  a  particle  travels  between  collisions.  Let  A  and  r  denote  the  mean-free  length  and  time, 
respectively,  characterizing  the  microscopic  particle  collisions.  Relevant  hydrodynamic  quantities 
are  the 

•  characteristic  flow  speed ,  u  ~ 

•  sound  speed,  cs  ~ 

•  shear  viscosity,  p  (and  the  kinematic  viscosity,  v  =  ^  ~  ^);  and, 

•  bulk  viscosity , 

The  relevant  dimensionless  quantities  are  the 

•  Knudsen  number ,  Kn,  defined  as  the  ratio  of  the  mean-free  path  to  the  characteristic  length 
scale  (Kn  =  £); 

•  Strouhal  number,  Sh,  defined  as  the  ratio  of  the  mean-free  time  to  the  characteristic  time  scale 
(Sh  =  f ); 

•  Mach  number,  M,  defined  as  the  ratio  of  the  characteristic  velocity  to  the  sound  speed  (M  =  — ); 

•  Reynolds  number,  Re,  defined  as  the  ratio  of  the  product  of  the  characteristic  velocity  times 
characteristic  length  to  the  kinematic  viscosity  (Re  =  ^  ~  Kn^’  and> 

•  fractional  mass  density  variation, 
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An  Analytical  Method  for 
Predicting  Macroscopic  Behavior 


Table  2.1:  Lattice-Gas  Model  Symbols 


Symbols 

Names 

i 

microscopic  cell  size 

4 

mesoscopic  cell  size 

T 

time  unit 

m 

single  particle  mass 

c 

velocity  unit  (-) 

D 

spatial  dimension 

B 

lattice  coordination  number 

a 

directional  index  (1,2,. . . ,  B) 

i,j,k,l 

spatial  indices 

& ai 

displacement  vectors 

2.1  Mesoscopic  Scale 

2.1.1  Derivation  of  the  Linearized  Kinetic  Equation 

Define  the  occupancy  probability  of  a  local  state,  fa  =  (na),  as  an  ensemble  average  over  inde¬ 
pendent  microscopic  realizations  of  the  lattice-gas  system.  In  the  mean-field  limit  approximate,  the 
particle  occupancies  /a’s  are  considered  uncorrelated  at  all  times.  Neglecting  the  particle-particle 
correlations  is  known  as  the  Boltzmann  molecular  chaos  assumption  (or  the  Stosszahlansatz) .  This 
approximation  allows  us  to  estimate  the  macroscopic  quantities,  such  as  the  damping  constants,  and 
the  shear  and  bulk  viscosities  (and  the  equation  of  state  in  multiphase  lattice  gases).  We  describe 
the  mesoscopic  dynamics  of  a  lattice-gas  system  at  the  mesoscopic  scale  by  the  lattice  Boltzmann 
equation 

fair  +  £ea,t  +  t)  =  fair, (2-1) 

where  the  momentum  index  a  ranges  from  0  to  B  —  1  specifying  the  local  states  at  each  site  of 
the  lattice.  The  superscript  on  the  collision  term,  fl”eso,  indicates  that  its  functional  form  is  to  be 
expressed  at  the  mesoscopic  scale.  The  asterisk  on  occupation  probability,  /*,  in  the  collision  term 
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Table  2.2:  Glossary  of  Single-Speed  Lattice  Gas  Variables  (L=l) 


Variables 

Names 

na 

microscopic  number  occupation  variable  for  a  local  state 

P 

microscopic  momentum 

microscopic  collision  operator 

s 

microscopic  configuration  of  particles  at  a  site 

TSS' 

transition  matrix 

fa 

mesoscopic  probability  of  occupation  of  a  local  state,  (na) 

Qmeso 

^  La 

mesoscopic  collision  operator 

Jab 

Jacobian  matrix, 

\a) 

eigenkets  of  J 

eigenvalues  of  J 

Cs 

macroscopic  sound  speed 

P 

macroscopic  pressure 

P 

macroscopic  density 

V 

macroscopic  velocity 

n  ij 

macroscopic  momentum  flux  density  tensor 

V 

shear  viscosity 

V 

kinematic  viscosity,  - 

indicates  that  fl“eso  depends  on  all  the  local  states  in  the  system  and  not  just  fa.  In  the  mean- field 
limit ,  f !“““(/*)  =  fl”f(/o,  fi,  ■  ■  ■ ,  fzs—i)  depends  only  on  the  local  states  at  the  lattice  site  r. 

The  parameters  t  and  r  are  the  cell  size  and  the  update  time,  respectively,  of  the  spacetime 
lattice  at  the  mesoscopic  scale.  The  parameter  B  equals  the  number  of  local  states  at  a  site  in  a 
Bravais  lattice.  The  vectors  ea  are  displacement  vectors  whose  length  corresponds  to  the  speed  of 
the  particles,  since  each  local  state  may  be  occupied  by  a  particle  with  momentum  p  =  mcea. 

2.1.2  The  Approach  to  Steady-State  Equilibrium 

The  system  is  said  to  be  in  steady-state  equilibrium  (which  may  also  be  called  thermodynamic 
equilibrium )  when  the  collision  term  in  the  lattice  Boltzmann  equation  vanishes 

nr°(/;q)  =  o,  (2.2) 

where  /°q  is  equilibrium  occupation  probabilities.  Therefore,  at  steady-state  equilibrium,  the  occu¬ 
pancy  probabilities  are  unchanging  over  time.  The  distribution  along  the  momentum  directions  of 
the  particle  occupancies  are  uniform,  so  the  local  configurations  are  perfectly  symmetric,  and  fl“eso 
cannot  cause  any  further  changes. 

Let  us  predict  the  non-equilibrium  behavior  of  the  lattice-gas  system  when  it  is  nearby  and 
approaching  steady-state  equilibrium.  Since  continuous  macroscopic  fields  for  the  mass  and  momen¬ 
tum  densites  are  defined  for  the  lattice-gas  system  in  the  continuum  limit,  by  Equations  (2.34)  and 
(2.35),  we  can  characterize  the  system  using  the  dimensionless  quantities  traditionally  used  to  char¬ 
acterize  fluid  systems.  Given  the  law  of  similarity1 ,  if  the  macroscopic  scale  behavior  of  the  lattice 
gas  is  fluid-like,  then  it  may  be  compared  to  a  natural  fluid  characterized  by  the  same  dimensionless 
quantities. 

Several  dimensionless  quantities  (the  Knudsen,  Strouhal,  Mach,  and  Reynolds  numbers,  and 
fractional  mass  density  variation)  allow  us  to  quantify  how  close  the  system  is  to  steady-state  equi¬ 
librium.  At  steady-state  equilibrium  all  the  dimensionless  numbers  vanish  (only  the  Mach  number 


1  See  page  56  of  Fluid  Mechanics  by  Laudau  and  Lifshitz  [96]. 


2.1.  MESOSCOPIC  SCALE 


13 


may  be  nonzero  at  equilibrium  if  there  is  a  global  uniform  background  flow,  however  this  can  be 
avoided  by  an  appropriate  choice  of  the  Galilean  frame-of-reference).  In  nonequilibrium  situations 
far  from  steady-state,  Kn,  Sh,  M,  and  ^  are  order  unity. 

Hydrodynamic  behavior  is  attained  in  the  long  wavelength  limit  where  Kn  and  Sh  are  close  to 
zero.  Viscous  hydrodynamic  behavior  is  attained  in  the  long  wavelength  limit  when  Sh  ~  Kn2 
and  —  ~  Kn.  This  is  called  diffusive  ordering  which  is  characteristic  of  random  walk  processes.2. 
Incompressible  viscous  hydrodynamics  occurs  when  we  also  have  M  ~  Kn  so  that  Re  ~  0(1)  and 
^  ~  Kn2.  A  procedure  for  linearizing  the  mesoscopic  Boltzmann  equation,  and  comparing  the 
resulting  dispersion  relations  to  the  solution  of  the  effective  field  theory  equations  of  motions,  is 
given  below  in  Section  2.1.3.  The  procedure  involves  a  series  expansions  in  The  standard  Mach 
number  expansion  of  the  probability  of  occupancy  is  used  [3].  Then,  a  Chapman-Enskog  procedure, 
given  in  Section  2.2.2,  which  is  necessary  for  the  derivation  of  the  macroscopic  equations  of  motion 
involves  perturbative  expansions  in  Kn  and  Sh,  given  in  Section  2.2.3. 

At  t  =  oo,  an  infinite  lattice-gas  system  completely  relaxes  to  steady-state  equilibrium,  where 
the  mass  density  field  is  uniformly  constant.  The  steady-state  equilibrium  occupation  probability, 
denoted  by  d ,  of  every  local  state  are  all  the  same  fa(x,oo)  =  d ,  for  all  a  and  all  x.  For  a  lattice 
of  finite  size,  the  number  of  phasespace  points  is  also  finite,  although  extremely  large.  The  number 
of  phasespace  points  equals  2BV ,  where  B  is  the  number  of  local  states  per  site  and  V  is  the  total 
number  of  sites.  Hence  the  Poincare  recurrence  time,  which  is  the  number  of  phasespace  points  of 
a  closed  loop  trajectory,  is  also  finite  and  the  state  of  the  finite-lattice  gas  system  is  not  defined  at 
t  =  oo.  Hence,  we  may  instead  say  that  the  lattice-system  has  completely  relaxed  to  a  steady  state 
on  a  time  scale  much  much  larger  than  the  characteristic  time  (gjj-)  for  the  largest  hydrodynamic 
scale  fluctuation. 

2.1.3  Linearized  Lattice  Boltzmann  Equation 

Our  first  step  towards  analyzing  the  nonsteady-state  behavior  of  the  system  will  be  to  expand 
jnoneq  ^  We  write  the  occupancy  probability  as  a  constant  part  (d  =  -*g)  and  a  fluctuating 

part  6 fa  (x ,  t)  <  d 

fa(x,t)  =  d  + 6fa(x,t).  (2.3) 

We  can  also  expand  the  collision  term  about  steady-state  equilibrium 

=  fi“ 60(/D  +  «Cao.  (2.4) 

Now  the  first  term  on  the  R.H.S.  vanishes  according  to  Equation  (2.2)  and  the  second  term  on  the 
R.H.S.  arises  because  of  fluctuations  in  the  probability  of  occupancies  of  all  the  local  states  in  the 
entire  system. 

At  the  mesoscopic  scale,  we  regard  each  of  the  occupancy  probabilities,  fa,  as  a  continuous 
variable.  The  basic  approach  is  that  Lla  is  a  continuous  and  differentiable  function  of  the  occu¬ 
pation  variables.  With  this  understanding,  using  the  chain  rule,  we  can  write  the  collision  term 
Equation  (2.4)  as  follows 


s(/:q)  =  =  J2 


b—1 


dfb 


Sfb  +  0(6f). 


f*=d 


(2.5) 


Using  Equations  (2.3)  and  (2.5),  we  can  write  the  lattice  Boltzmann  equation  Equation  (2.1)  in 
linearized  form 

Sfa(x  +  £sea,t  +  T)  =  S  fa  (x ,  t)  +  JabSfb(rff),  (2.6) 

2  See  Figure  3.9  using  the  classical  lD3Px  model  to  illustrate  diffusive  ordering  that  is  characteristic  of  lattice-gas 
fluids.  Furthermore,  a  tagged  particle  in  a  lattice  gas  undergoes  a  random  walk,  and  the  observed  diffusive  behavior 
of  tagged-particles  in  lattice-gas  simulations  agrees  well  with  analytical  results  [97,  29] 
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where  the  Jacobian  of  the  collision  term  is  defined  as 


Jab  — 


dfb 


L=d 


(2.7) 


Let  fa(k,  to)  denote  the  discrete  Fourier  transform  of  the  occupation  probability  /a(cc,  t).  Then  taking 
the  discrete  Fourier  transform  of  the  linearized  Boltzmann  equation  Equation  (2.6)  we  obtain  the 
following  characteristic  equation 

ei(wS+“T)5/a(£,w)  =  Sfa{k,w)  +  JabSfb(k,  to),  (2.8) 


which  we  rewrite  as 

Sfb(k,  w)  =  0.  (2.9) 

Therefore,  we  have  the  following  matrix  equation 


fei(esea.k+uT)  _  Sab  _  Jab 


M6f  =  0,  (2.10) 

where  Sf  =  (Sfo,  8f\, . . . ,  Sfs—i)  and  the  components  of  the  square  matrix  M  are 

Mab  ee  ^  Sab  _  Jab  (2.H) 

Solving  Equation  (2.10)  gives  us  the  dispersion  relations  for  the  system  obeying  what  is  called 
generalized  hydrodynamics.  The  generalized  hydrodynamics  for  classical  lattice-gas  systems  have 
been  previously  worked  out  [30,  98].  The  development  given  here  for  the  lattice-gas  system  follows 
Das  and  Ernst’s  treatment  of  a  classical  lattice-gas  system  [30].  However,  in  the  present  treatment,  I 
do  not  use  differential  point  form  notation  for  mesoscopic  fields  (since  technically  this  is  unwarranted 
and  allowed  only  in  the  continuum  limit).  Instead,  to  be  absolutely  rigorous,  I  have  applied  the 
discrete  Fourier  transform  to  the  mesoscopic  field  to  obtain  Equations  (2.10)  and  (2.11).  So  up  to 
this  point  in  the  analytical  treatment,  I  have  not  invoked  the  continuum  limit. 


2.1.4  Dispersion  Relations 

To  solve  Equation  (2.10)  for  the  dispersion  relation  to  =  u(k),  we  must  find  the  B  roots  of  the 
secular  determinant  of  the  matrix  M  [30,  98].  The  solution  can  be  found  analytically  for  the  lD3Px 
lattice  gas  and  numerically  for  more  complicated  lattice  gases.  In  general,  there  are  two  types  of 
long  wavelength  excitations  {k  — >  0),  and  they  are  called  hard  kinetic  modes  and  soft  hydrodynamic 
modes. 

In  the  long-wavelength  limit,  the  kinetic  modes  are  nonvanishing  at  first  order.  The  kinetic 
modes  decay  rapidly  in  the  lattice-gas  system  because  of  a  positive  imaginary  part  in  the  eigenvalue 
spectrum  of  u>,  (Im(w)  >  0),  at  k  =  0.  In  contrast,  the  soft  hydrodynamic  modes  decay  over  a 
long-time  scale.  They  are  associated  with  eigenvalues  that  vanish  in  the  long-wavelength  limit, 
(Im(u;)  =  0),  at  k  =  0.  These  vanishing  eigenvalues  in  turn  are  associated  with  the  conserved 
quantities  in  the  lattice-gas  system. 

In  the  long- wavelength  limit,  k  ~  0,  the  dispersion  relation  for  the  lD3Px  lattice-gas  corresponds 
to  a  damped  sound  mode 

ui(k)  =  ±csk  +  iT(p)k2.  (2.12) 

The  real  part  of  ui  is  linear  in  the  wave  number.  Since  Re(w)  is  linear  in  the  wave  number,  the  sound 
mode  excitation  propagates  at  the  sound  velocity  corresponding  to  the  slope,  which  is  denoted  here 
by  cs.  Furthermore,  the  sound  mode  is  damped  in  the  viscous  hydrodynamic  regime  characterized 
by  diffusive  ordering  where  the  dispersion  relation  is  parabolic  in  wavenumber,  Im(u;)  ~  0[k2). 
In  general,  in  a  single-speed  lattice  gas,  the  decay  of  the  hydrodynamic  modes  depends  on  shear 
viscosity  and  sound  damping  (there  is  no  mode  related  to  bulk  viscosity). 
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2.1.5  Criterion  for  Deviations  from  Local  Equilibria 


At  local  equilibrium,  we  assume  that  the  occupancies  of  the  local  states  at  each  site  in  the  system 
are  isotropic.  This  is  called  the  subsonic  limit.  A  stronger  definition  of  the  subsonic  limit  is  when 
the  fractional  variation  of  the  occupancy  probabilities,  \6fa\/ at  all  lattice  sites,  are  assumed  to 
be  uniformly  distributed  along  all  the  momentum  directions,  1  <  a  <  B.  Hence,  the  criterion  for  a 
fractional  mass  density  variation  on  the  order  of  the  Knudsen  number,  can  be  expressed  as3 


\Sfa\ 

f?  ~  BL' 


(2.13) 


The  maximum  size  of  the  simulation  volume  is  limited  by  the  amount  of  physical  computational 
resources  available.  In  a  microscopic  lattice-gas  simulation,  the  occupation  probabilities,  fa ,  must 
be  determined  by  either  partitioning  the  maximum  size  simulation  into  ensemble  realizations  or 
coarse-grain  blocks.  Let  us,  for  the  moment,  revisit  the  subject  of  averaging  over  the  microscopic 
quantities  to  obtain  the  mesoscopic  values.  There  are  more  details  to  discuss. 

In  ensemble  averaging,  many  realizations  of  the  lattice-gas  system,  which  are  identical  at  the 
macroscopic  scale,  are  computed  independently.  Measurements  are  separately  made  from  each  re¬ 
alization,  and  all  the  resulting  measurements  are  then  averaged.  For  example,  the  state  of  the  ath 
qubit  is  measured  for  each  copy  of  the  system,  which  results  in  a  series  of  l’s  and  0’s,  and  the  average 
value  is  an  estimate  of  fa.  In  coarse-grain  averaging,  one  measures  the  occupancy  of  all  the  local 
states,  occupied  by  a  particle  with  momentum  mcea,  at  all  the  sites  within  a  spacetime  block  of  a 
large  microscopic  system.  Again,  this  results  in  a  series  of  l’s  and  0’s,  which  are  then  averaged  to 
estimate  fa(X,T),  where  say  X  and  T  denote  the  coordinates  of  the  centroid  of  a  spacetime  block 
within  the  superlattice. 

In  either  case,  whether  ensemble  or  coarse-grain  averaging  is  employed,  all  the  available  com¬ 
putational  resources  are  expended.  However,  the  numerical  results  can  be  quite  different,  for  two 
reasons:  because  lattice-gas  systems  obey  diffusive  ordering  and  because  of  renormalization  effects 
arising  from  particle-particle  correlations. 

Let  us  first  consider  the  consequences  of  diffusive  ordering.  If  one  doubles  the  system  size, 
L  — >  2 L,  one  must  quadruple  the  simulation  time,  T  — >  AT,  to  evolve  to  a  macroscopic  state 
similar  to  the  one  obtained  by  running  a  simulation  of  size  L  for  time  T.  Consequently,  if  the 
fixed  amount  of  computational  resource  is  partitioned  to  do  ensemble  averaging,  then  many  “small” 
systems  are  simulated  which  “rapidly”  relax  towards  steady-state  equilibrium.  If  the  fixed  amount 
of  computation  resource  is  partitioned  to  do  coarse-grain  averaging,  then  one  “large”  system  is 
simulated  which  “slowly”  relaxes  towards  steady  state.  Therefore,  estimates  can  be  made  more 
quickly  using  ensemble  averaging,  but  only  in  those  situations  were  particle-particle  correlations  can 
be  neglected. 

This  brings  us  to  the  next  issue  of  renormalization.  In  a  large  system  simulation,  there  is 
sufficient  time  for  many  particle  collision  events  to  occur  allowing  the  particle  occupancies  to  become 
correlated.  These  particle-particle  correlations,  in  certain  situations,  may  have  an  appreciable  effect 
on  the  value  of  the  transport  coefficients  [88].  One  enumerates  all  connected  diagrams  corresponding 
to  the  pathways  by  which  outgoing  particles,  initially  correlated  by  a  collision,  move  through  the 
system,  interact  with  other  particles,  and  eventually  return  as  incoming  particles  to  a  final  collision 
event.  Each  connected  diagram  corresponds  to  a  term  in  an  asymptotic  series  expansion  of  the 
collision  operator,  which  is  summed  to  give  a  renormalized  collision  operator.  If  the  ultimate  aim  of 
is  to  estimate  a  transport  coefficient,  which  strongly  depends  on  particle-particle  correlations,  then 
coarse-graining  averaging  must  be  used  [99,  88]. 

3  Suppose,  as  a  concrete  example,  that  an  initial  nonequilibrium  state  of  the  system  is  chosen  with  a  characteristic 
feature  size  on  the  order  of  say  one  hundred  lattice  grid  units,  L  ~  1 00/:.  Next,  suppose  the  mean-free  path  length 
is  on  the  order  of  the  size  of  a  single  primitive  lattice  cell,  A  ~  t.  If  the  fractional  mass  density  variation,  which 
must  be  on  the  order  of  the  Knudsen  number  is  ^  ~  L  ^  0.01,  then  the  lattice  gas  would  accurately  model  the 
dynamical  fluid  behavior  in  the  regime  of  incompressible  viscous  hydrodynamics.  Continuing  the  example,  if  B  =  6 
then  |<5/„|//'q  ~  0.002. 
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In  either  case,  it  is  necessary  to  satisfy  the  criterion  that  ( Sfa/fa q)  ~  0( Kn)  or  smaller  for  all  a. 
If  this  criterion  is  satisfied,  the  linearized  lattice  Boltzmann  equation  Equation  (2.6)  can  accurately 
describe  Navier-Stokes  hydrodynamics.  And  this  is  a  stronger  requirement  then  ^  ~  0( Kn).  The 
requirement  that  ( Sfa/fa q)  ~  0(Kn)  for  viscous  hydrodynamics  (or  the  more  stringent  requirement 
that  (< 5/a//aq)  ~  0( Kn2)  for  incompressible  viscous  hydrodynamics)  implies  a  lower  bound  for  the 
number  of  states  used  in  an  ensemble  average  or  for  the  minimum  size  of  the  spacetime  block  used  in 
a  coarse-grain  average.  All  these  considerations,  usually  applied  to  the  classical  lattice-gas  method, 
are  also  relevant  to  quantum  lattice-gas  simulations  [100,  101],  which  are  not  reviewed  here. 


2.2  Macroscopic  Scale 

2.2.1  Eigensystem  of  the  Linearized  Collision  Operator 

In  the  long  wavelength  (k  — »•  0)  limit,  the  characteristic  equation  Equation  (2.9)  reduces  to  the 
simple  form 

[(e“T  —  1)1  —  J]  <Sf  =  0.  (2.14) 

Expanding  to  first  order  in  Sh  (second  order  in  e)  this  become  the  eigenvalue  equation 

J  Sf  =  ujrSf  +  0(£3).  (2.15) 

Therefore,  in  the  long- wavelength  and  low-frequency  limits,  the  eigenvalues  of  J  determine  possible 
values  for  u>  and  in  turn  the  hydrodynamic  and  kinetic  behavior  of  the  lattice-gas  system.  This 
eigenvalue  problem  is  analytically  solvable,  without  the  need  for  any  numerical  treatment  as  is  needed 
for  finding  the  fc-dependent  roots  of  the  secular  determinant  of  the  matrix  M  in  Equation  (2.11). 
Consider  the  following  eigenvalue  equation 


Jabtf  = 


(2.16) 


with  eigenvectors  and  eigenvalues  na,  where  a  =  1, . . . ,  B.A  Let  us  see  why  there  will  be  as  many 
zero  eigenvalues  as  there  are  conserved  quantities  in  the  lattice  gas  dynamics.  For  convenience,  we 
will  use  ket  notation  where  jet)  =  (£f ,  . . . ,  £g)  and  | Sf)  =  {Sfi,  6/2, . . . ,  Sfs)-  We  can  write  J  as 

follows 

B 

j  =  hh  (2-19) 

Q!=l 

so  Equation  (2.15)  becomes 

B 

Y.  Ka\a){a\Sf)  =  ujr\6f).  (2.20) 

a=l 


All  the  scalars  (a\Sf<'1'))  for  which  na  =  0  have  no  effect  on  the  dynamics  since  J\Sf)  =  0  and 
so  correspond  to  the  conserved  quantities  of  the  system.  The  set  of  eigenvectors  with  degenerate 
eigenvalue  of  zero  span  what  is  called  the  hydrodynamic  space ,  which  I  denote  by  H.  The  remaining 
set  of  eigenvectors  (with  nonzero  eigenvalues)  span  what  is  called  the  kinetic  space ,  which  I  denote 

4  The  problem  is  simplified  if  J  is  circulant  [  ].  The  components  of  J  can  be  specified  by  the  difference  of  the 
indices,  J,,},  =  Ja—b ■  Hence,  we  make  the  ansatz  that  the  eigenvectors  £“  have  the  following  form 

^=e2™aa/B.  (2.17) 

Then  inserting  Equation  (2.17)  into  Equation  (2.16)  and  taking  m  =  a  —  6,  gives  a  solution  for  the  eigenvalues 


m= 1 


,2'Kimcx.  /  B 


(2.18) 
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as  1C.  Therefore  J  can  be  explicitly  written  as  a  linear  combination  over  eigenvectors  in  the  kinetic 
space 

J=J2K*  hh-  (2-21) 

aelC 

In  an  athermal  system,  there  are  1  +  D  conserved  quantities,  the  mass  plus  the  momentum  for 
each  dimension  of  the  space.  So  LL  is  a  D  +  1  dimensional  space  and  /CisaB  —  D  ^  1  dimensional 
space.  Let  us  denote  the  kinetic  eigenkets  as  follows,  \D  +  2),\D  +  3), . . .  ,\B),  which  span  the  kinetic 
subspace  1C.  Then,  the  generalized  inverse  of  J  is  defined  over  K.  as  follows 


J1  =  (\D  +  2)  \D  +  3)---\B)) 

" - V - ' 

BxK  matrix 


/  0 
r  «£>+2 

0  1 

KD+ 3 

o  o 

/  (D  +  2  \ 
(D  +  3 

VO 

V  (B\  J 

KxK  matrix 


KxB  matrix 


(2.22) 


From  Equation  (2.22),  it  follows  by  construction  that 

J_1|a)  =  —  |  a), 


(2.23) 


for  |a)  £  1C. 

Within  /C  there  exists  the  viscous  subspace ,  V  C  1C  characterized  by  the  degenerate  eigenvalue  nv . 
Because  the  collisional  process  is  invariant  under  the  finite  point-group  symmetries  of  the  Bravais 
lattice,  consequently  there  is  a  lack  of  preference  in  direction  for  momentum  diffusion  in  the  system. 
Hence,  there  is  a  subspace  of  K.  characterized  a  degenerate  eigenvalue  which  contributes  positively 
to  the  shear  and  bulk  viscosities. 

A  ket  |e,;ej)  may  be  formed  from  the  dyadic  product  eaieaj.  The  ket  |  CjCj)  resides  in  V.  It  is  an 
eigenket  of  the  generalized  inverse  of  the  Jacobian  of  the  collision  operator,  J-1,  and  has  eigenvalue 
-A.  That  is,  for  a  hydrodynamic  lattice-gas  fluid 

J~x\eiej)  =  — |  eiCj).  (2.24) 

The  identity  Equation  (2.24)  will  be  needed  in  the  following  section. 


2.2.2  Chapman-Enskog  Expansion 


The  characteristic  equation  Equation  (2.9)  of  the  linearized  lattice  Boltzmann  equation 


sea-k-\-u)T 


—  l)  Sab 


Jab 


Sfb(k,u)  =  0 


is  an  approximate  description  of  the  mesoscopic  particle  dynamics  since  the  collision  term  on  the 
R.H.S.  has  been  expanded  to  first  order  about  the  equilibrium  value  of  the  occupation  probability. 
In  this  approximation,  the  collision  term,  Jab,  acts  as  a  linear  operator  on  the  local  configuration  <5f . 
I  would  now  like  to  expand  the  L.H.S.  of  this  characteristic  equation.  To  do  so,  let  us  use  e  as  a  small 
expansion  parameter,  e  <C  1.  In  the  viscous  hydrodynamic  regime,  this  expansion  parameter  is  the 
Knudsen  number,  Kn  ~  ea  •  k  =  £s\k\  ~  e.  And,  because  of  diffusive  ordering,  the  Strouhal  number 
is  Sh  ~ur~  e2.  We  expect  Equation  (2.9)  is  an  appropriate  description  of  the  mesoscopic  dynamics 
so  long  as  the  nonequilibrium  occupation  probabilities  are  close  enough  to  their  equilibrium  values 
so  that  the  action  of  the  linearized  collision  term,  Jab ,  is  sufficient  to  cause  any  such  nonequilibrium 
configuration  to  relax  back  to  an  equilibrium  configuration. 

Begin  by  expanding  Equation  (2.9)  to  first  order  in  e 


( i£sea  *  kSab  Jab)Sfb  —  0- 


(2.25) 
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The  basic  approach  is  that  the  deviations,  Sfa ,  of  the  occupation  probability  can  be  expanded  in 
powers  of  e 

5fa  =  Sf^  +Sf^  +0(e2),  (2.26) 

where  Sfa0'1  ~  e  and  Sf ^  ~  £2.  The  superscript  on  Sf^  denotes  that  it  is  a  deviation  from  the 
steady-state  equilibrium  d  due  to  bulk  motion  of  the  fluid.  The  superscript  on  Sf^  denotes  that  it  is 
a  deviation  due  to  spatial  gradients  in  the  bulk  profile.  We  require  that  the  largest  Mach  number  at 
any  point  in  the  system  must  be  small,  M  <C  1,  that  Sf^  ~  M  ~  e.  We  can  insert  Equation  (2.26) 
into  the  first  order  e-expansion  of  the  characteristic  equation  Equation  (2.25),  JabSfb°^  =  0,  and  we 
equate  the  two  0{e2)  terms 

Usea-kSf^  =  Jab8f^.  (2.27) 

Since  J  has  a  well-defined  generalized  inverse,  we  can  invert  the  Jacobian  matrix  according  to 
Equation  (2.22)  to  solve  for  the  second  order  correction  to  the  occupation  probability 

<S/(1)  =iU~b1eb-k6^0).  (2.28) 

Therefore,  using  the  basic  approach  of  the  Chapman-Enskog  expansion  Equation  (2.26),  we  have 
the  result  that 

Sfa  =  [Sab  +  iU~bh  ■  k\S /b(0)  +  0(e3).  (2.29) 

In  the  continuum  limit,  we  are  justified  in  taking  the  inverse  Fourier  transform  of  Equation  (2.29), 
which  gives  the  fluctuating  part  of  the  non-equilibrium  probability  occupancy  in  differential  point 
form 

6fa(x,t)  =  [<U  +  U~bh  •  V]<5/b(0)(x,t)  +  0(e3).  (2.30) 

Then  using  Equation  (2.3),  the  probability  of  occupancy  in  the  continuum  limit  is 

fa(x,t)  =  [Sab  +  ■  V]/6(0)(5?,t)  +  (D(e3).  (2.31) 

We  can  insert  the  standard  Mach  number  expansion  [3]  into  Equation  (2.31).  After  some  algebraic 
manipulation,  the  result  is 

fa  =  d[  1  +  ^eaiVi  +  QaijViVj  +  TDJ~b1ebiebjdzvj}  +  0{M3).  (2.32) 

Using  the  identity  Equation  (2.24),  that  J~^ebiebj  =  eaieaj ,  Equation  (2.32)  becomes 

fa  =  d[  1  +  —eaiVi  +  +  QaijViVj  -  —eaieajdiVj]  +  0(M3).  (2.33) 

C  ZiC  txrq 

The  approximation  Equation  (2.33)  is  a  good  one  provided  several  conditions  are  met: 

1.  the  ratio  of  the  superlattice  cell  size  to  the  characteristic  scale  length  of  the  small  hydrody¬ 
namics  fluctuation  is  close  to  zero,  ^  ~  0  (satisfied  in  the  continuum  limit) 

2.  the  angular  distribution  of  particles  along  momentum  directions  is  close  to  an  isotropic  one 

3.  the  flow  is  subsonic  (M  -C  1) 

4.  spatial  gradients  are  small  (Kn  ~  0  and  ^  is  small) . 
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2.2.3  Derivation  of  the  Continuum  Equations 

The  local  mass  density  and  the  momentum  density  at  x  and  t  can  be  expressed  in  terms  of  the 
occupancy  probability,  fa(x,t),  is 


B 

p{x,  t)  =  lim  /a(ir,f)  (2.34) 

•ts— >■() 

a=  1 
B 

p(x,t)vi(x,t)  =  lim  ^2  me2 eai  fa{x,t). 

“  a=l 

(2.35) 


The  mass  and  momentum  densities  are  considered  “macroscopic”  field  quantities.  They  are  only  well 
defined  in  the  continuum  limit ,  where  the  primitive  cell  size  of  the  lattice  approaches  zero.  However, 
for  practical  considerations,  they  are  approximated  by  high  resolution  grids  with  small  but  finite 
cell  size. 

The  derivation  of  the  continuum  equations  of  motion  at  the  macroscopic  scale  is  carried  out  in 
this  section.  The  method  of  derivation  is  outlined  by  these  following  few  steps: 

1.  Expand  the  mesoscopic  Boltzmann  equation  to  first  order  in  time  and  second  order  in  space.5 

2.  In  the  continuum  limit,  calculate  the  first  and  second  moments  of  the  mesoscopic  Boltzmann 
equation  .6 

3.  Insert  the  mesoscopic  occupancy  probability  given  Equation  (2.33)  into  the  moment  equations 
obtained  in  Step  2. 


After  some  algebraic  manipulations,  we  thereby  obtain  an  approximation  of  the  equations  of  motion 
that  serve  as  an  effective  field  theory  at  the  macroscopic  scale.  Because  of  diffusive  ordering,  the 
result  is  that  the  macroscopic  equations  of  motion  are  a  set  of  coupled  parabolic  partial  differential 
equations.  In  the  present  derivation,  I  do  not  give  a  multi-scale  analysis  such  as  the  one  carried 
out  by  Frisch  et  al.  [3]  in  their  treatment  of  2  and  3  dimensional  lattice-gas  hydrodynamics.  This 
omission  is  justified  because  only  a  single  time-scale  is  needed  for  most  lattice-gas  systems  since 
the  transport  coefficients  are  very  large.  That  is,  there  is  little  or  no  separation  between  the  short 
time  scale  associated  with  sound  mode  excitations  and  the  longer  time  scale  associated  with  viscous 
mode  excitations  arising  from  momentum  diffusion.  Viscous  damping  in  lattice-gas  fluids  is  observed 
over  relatively  short-time  scales  and  therefore  significantly  affects,  and  is  mixed  in  with,  sound  wave 
propagation  and  convection.  Nevertheless,  at  the  end  of  this  section,  I  will  divide  the  effective  field 
theory  into  two  sets  of  equations  that  apply  at  short  and  long  time  scales. 

We  determine  the  macroscopic  equations  of  motion  using  the  lattice  Boltzmann  equation 


dfa 

dt 


lim 

T->0 


ft" 


(2.36) 


In  consideration  of  diffusive  ordering,  we  expand  the  L.H.S.  of  this  equation  to  first-order  in  time 
and  second-order  in  space 

F)  f  P ^  Qmeso 

-fey  +  cea  •  V/a  +  -^-(eQ  •  V)2/a  +  0(Kn3,Sh2)  =  lim  .  (2.37) 

ot  zt  t->o  r 

la->  0 

5  Only  a  first  order  time  derivative  is  needed  because  of  the  long-time  scales  associated  with  viscous  damping. 
Time  and  spatial  scales  are  related  parabolically  (T  ~  L 2  or  e2  ~  Sx2  ~  St)  in  lattice-gas  systems. 

6  This  is  done  because  for  each  additive  conserved  quantity  of  the  dynamics,  a  macroscopic  field  is  expressed  as  a 
moment  of  the  mesoscopic  field  of  probability  of  occupancies.  In  the  present  case,  p  and  v  are  expressed  in  terms  of 
the  fa  s  according  Equations  (2.34)  and  (2.35). 
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This  is  called  the  lattice- Boltzmann  equation ,  where  12“eao  is  defined  by  Equation  (2.5).  Since 
Y^a  fi”eso  =  0,  the  zeroth  moment  of  Equation  (2.37)  is 


dt  faJ  +  d i  |  rnc^eaifa  j  +  T^didj  (  eaieajfa  J  +  £>(Kn3, Sh2)  =  0. 
Using  the  identities  [1] 

£< 

a=  1 

B 

'y  ( eaea 

a=l 
B 

y  ( eaeaea 

a= 1 
B 

£eae 


e„  =  0 


(2.38) 


-A<2) 

D 


=  0 


B 


£>(£>  +  2) 


A<4> 


(2.39) 

(2.40) 

(2.41) 

(2.42) 


the  corrected  occupancy  probability  Equation  (2.33),  along  with  definitions  for  the  mass  density 
Equation  (2.34)  and  momentum  density  Equation  (2.35),  this  reduces  to  a  mass  continuity  equation 
in  the  long- wavelength,  low-frequency,  and  subsonic  limits 


dtp+diipvi)  +C>(Kn3,Sh2,M3)  =  0. 
Since  )T)n  eaiLl™so  =  0  too,  the  first  moment  of  Equation  (2.37)  is 


(2.43) 


dt  |  uic  y  ^  eaifa  J  +  ()j  |  me  y  ( ea i.eajfa  J  +  ^ ^_djd^  ^  ^  eai€ajCakfa J  +  £>( Kn  ,  Sh  )  0. 

(2.44) 


This  reduces  to  Euler’s  equation  in  the  long- wavelength,  low-frequency,  and  subsonic  limits 

dt(pvi)  +  djHij  +  0(Kn3,  Sh2,  M3)  =  0,  (2.45) 


where  the  momentum  flux  density  is 

II ij  —  Pij  H-  gpViVj 


(A  fl  ,  Ps  P  a 
(£>  +  2)TKvdj Vi  +  2t  (£>  +  2)djVi 


(2.46) 


or 


n,;,  =  Pij  +  gpViVj  - 


pi] 


(£>  +  2 )t  kv 


(2.47) 


The  shift  of  is  a  constant  negative  contribution  to  the  shear  viscosity  by  a  lattice  effect.  With 
sound  speed  cs  =  identify  the  isotropic  pressure  tensor  as 


Pij  =  PCS  (  1  5  c2  £2  )  • 


(2.48) 


Finally,  inserting  Equation  (2.47)  into  Euler’s  equation  Equation  (2.45),  the  momentum  equation 
for  viscous  flow  is 


dt(pvi)  +  dj(gpViVj)  +  0( Kn3,  Sh2,M3)  =  -diP  +  r]d2Vi  +  (c  +  didjVj, 


(2.49) 
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with  shear  viscosity 

1  T  D  +  2\kv  2)  ’ 

and  bulk  viscosity 

pel  2D -I  f  1  1\ 

^  r  £>(£>  + 2)  ^  2/ 


(2.50) 


(2.51) 


With  small  Knudsen,  Strouhal,  and  Mach  numbers,  the  momentum  equation  Equation  (2.49)  ap¬ 
proximates  the  Navier-Stokes  equation  except  that  there  is  an  extra  density-dependent  factor,  g{d), 
in  the  convective  term  in  Equation  (2.49).  If  g  is  not  unity,  Galilean  invariance  is  destroyed.  It  is 
possible  to  alter  the  collision  term  in  the  lattice  Boltzmann  equation  so  that  g  =  1.  This  was  first 
done  by  H.  Chen,  S.  Chen,  and  Mattaeus  [59],  but  their  approach  violated  detailed  balance,  and 
consequently  the  algorthm  is  subject  to  unphysical  numerical  instabilities.  Unconditionally  numer¬ 
ically  stable  lattice-gas  methods  that  obey  the  principle  of  detailed  balance  have  also  been  found  to 
set  g  =  1  [5,  39]. 

For  lattice-gas  systems  with  low  viscosity  [102],  it  is  appropriate  to  consider  two  separate  ef¬ 
fective  field  theories  for  short-time  and  long-time  hydrodynamic  behavior.  The  short  and  long 
hydrodynamic  time  scales  are  denoted  by  t\  and  £2 ,  respectively.  A  multi-scale  formalism  is  used 
where  <9t  — >  <9ti  +  dt2  [3].  Then  the  effective  field  theory  defined  by  the  continuity  equation  Equa¬ 
tion  (2.43)  and  the  Navier-Stokes  equation  Equation  (2.45)  reduces  at  the  short-time  hydrodynamic 
scale  to  the  following  set 


d^p  +  ditpvi)  =  0  (2.52) 

dtl{pVl)  =  -djP.  (2.53) 

This  set  models  sound  wave  propagation  induced  by  pressure  gradients  in  a  compressible  fluid.  At 
the  long-time  hydrodynamics  scale,  Equations  (2.43)  and  (2.45)  reduce  to 

dt2p  =  0 

dt2(pi’i)  +  dj(gpviVj)  =  gd2Vi  +  (C  +  ^)  fydjVj. 


This  set  of  equations  describes  viscous  damping  and  convection  for  a  compressible  fluid. 


(2.54) 

(2.55) 
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Chapter  3 

A  One-Dimensional  Model  with 
Conserved  Mass  and  Momentum 


3.1  Microscopic  Rules 

Let  us  consider  the  simplest  lattice-gas  model  of  a  one-dimensional  fluid  with  two  conserved 
quantities.  This  is  called  the  IDSPx  lattice-gas  model.  This  model  was  first  studied  by  Qian  in  1990 
[103].  The  lattice  gas  is  one  dimensional  and  has  three  bits  per  site,  a  rest  particle  with  mass  two 
and  speed  ±1  particles  with  mass  one.  The  mass  and  momentum  at  a  lattice  site  is 

m  =  2uq  +  n±  +  ri2  and  px  =  n\  —  «2-  (3-1) 


m  =  2 

head-on 

rest 

Figure  3.1:  Head-on  collision  in  the  lD3Px  lattice-gas  model.  The  single  equivalence  class  has  m  =  2  and  px  =  0. 

There  are  two  local  configurations  both  with  m  =  2  and  px  =  0:  (1)  {no,ni,ri2}  =  {1,0,0}  and 
(2)  {no,  n,i,n,2}  =  {0, 1, 1}.  These  two  configurations  are  members  of  the  only  collision  set  (which  is 
called  an  equivalence  class).  An  equivalence  class  has  two  or  more  members.  Figure  3.1  illustrates 
the  equivalence  class  of  the  lD3Px  model.  Its  two  elements  are  the  configuration  of  two  head-on 
particles  {011}  and  the  configuration  with  a  single  rest  particle  {100}. 

Because  the  total  number  of  particles  and  the  total  momentum  must  be  conserved,  the  collision 
part  of  the  dynamics  can  only  permute  the  local  configurations.  The  collision  equation,  which  is 
applied  homogeneously  across  the  lattice,  can  be  expressed  in  terms  of  a  mapping  function,  U,  as 
follows 

s'  =  U(s),  (3.2) 

where  U  maps  2B  configurations  to  2s  new  configurations.  For  the  simple  lD3Px  lattice,  U  is 

C/({011})  =  {100} 

U{{  100})  =  {011}. 

If  a  configuration  s  is  not  a  member  of  an  equivalence  class,  then  U(s)  =  s.  In  other  words,  if  the 
incoming  state  is  not  a  member  of  an  equivalence  class,  then  the  outgoing  state  is  set  equal  to  the 
incoming  state.  To  speed  up  a  lattice-gas  simulation,  the  mapping  function,  U ,  may  be  precomputed 
before  the  simulation  and  accessed  in  lookup  table  fashion  during  the  simulation. 
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In  a  computer  implementation,  it  is  convenient  to  use  two  arrays  to  simultaneously  store  the 
states  s  and  s'.  Therefore,  in  Equation  (3.2),  data  in  the  array  that  stores  the  “incoming”  state,  s, 
is  transformed  by  the  action  of  the  lookup  table  U  (which  is  applied  homogeneously  over  the  entire 
array)  and  the  output  is  written  into  the  next  array  to  store  the  new  “outgoing”  state,  s' . 

It  is  conventional  to  write  the  collision  rule  in  terms  of  the  occupation  variables,  na  =  1  or  0, 
which  are  Boolean  values.  The  collision  rule,  expressed  for  an  individual  local  state,  is  written 

n’a(x,t)  =  na(x,t)  +  (3.3) 

where  the  collision  term  fla(n»)  =  ±1  or  0.  Writing  f la,(n*)  with  an  asterisk  subscript  on  n*  denotes 
that  the  collision  term  for  the  ath  local  state  depends  on  all  the  on-site  local  states.  It  is  conventional 
to  write  the  streaming  rule  in  terms  of  na  also 

na(x  +  £ea,t  +  T)  =  n'a(x,t).  (3.4) 

Combining  Equations  (3.3)  and  (3.4),  the  microscopic  transport  equation  is  therefore 

na{x  +  £ea,  t  +  r)  =  na{x,  t)  +  f la(n*).  (3.5) 

For  the  lD3Px  model,  the  lattice  vectors  are  eo  =  0,  ei  =  x,  and  e2  =  —x,  and  the  collision  term  is 
specified  by  the  single  function 


fi  =  rii  n2  (1  -  no)  —  n0  (1  -  ni)(l  -  n2).  (3.6) 

where  Qq  =  fh  and  =  —  Cl.  Then  explicitly  for  the  lD3Px  model,  the  microscopic  transport 
equation  Equation  (3.5)  is 

no(x,t  +  r)  =  no(x,  t)  +  Ll(x,  t)  (3-7) 

ni}2(x  ±  £,  t  +  t)  =  nit2(x,t)  —  Ll(x,t). 

Successive  snap  shots  of  the  microscopic  time  evolution  of  a  small  classical  lattice-gas  system  on 
a  one-dimensional  linear  grid  with  1024  sites  is  shown  in  Figure  3.2.  The  particles  (which  look  like 
lines  on  a  barcode)  are  rendered  in  black  on  the  right  side  of  the  diagram.  Each  cell  comprises  3  local 
states  and  hence  each  cell  can  hold  up  to  3  particles.  Bit  0  encodes  the  m  =  2  “rest”  particle  state, 
and  bits  1  and  2  encodes  states  with  m  =  1  particles  moving  with  speed  +1  and  -1,  respectively. 
Initially,  at  t  =  0,  most  of  the  particles  reside  on  the  left  half  of  the  lattice  and  successive  snapshots 
are  shown  at  every  128  time  steps.  A  spatial  coarse-grain  average,  with  block  size  of  128,  is  used  to 
determine  the  “mass  density  field”;  the  total  mass  per  site  is  plotted  on  the  left  side  of  the  figure. 
The  initial  state,  rendered  in  red,  is  a  sinusoidal  perturbation  with  a  background  density  of  p  =  \ 
and  an  amplitude  of  dp  =  ) .  (The  magnitude  of  5p  is  exaggerated  here  for  pedagogical  purposes;  it 
is  usually  a  few  percent  of  the  background  density.)  At  t  =  640,  the  situation  is  reversed  as  most 
of  the  particles  have  moved  to  local  states  on  the  right  half  of  the  lattice.  The  final  mass  density 
waveform  is  rendered  in  blue.  Each  microscopic  configuration  of  particles  rendered  in  this  figure 
corresponds  to  a  unique  state  of  the  lattice-gas  system.  Even  though  this  is  a  small  system  with 
only  1024  lattice  sites,  the  total  number  of  unique  states  available  is  still  an  extremely  large  number, 
23072.  The  sound  speed  is  clearly  less  than  unity,  because  for  the  wave  to  move  halfway  cross  the 
grid,  512  lattice  units,  it  takes  over  640  time  steps. 


3.2  Model  Analysis 

A  lattice  Boltzmann  equation  describes  the  dynamics  of  the  lD3Px  lattice-gas  system  at  the 
mesoscopic  scale.  The  mesoscopic  average  of  the  occupation  variable,  na(x,t),  is  the  probability  of 
occupancy 


fa(x,t )  =  (na{x,t)}. 


(3.8) 
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Figure  3.2:  Successive  snapshots  of  the  microscopic  state  of  a  classical  lD3Px  lattice-gas  system  with  a  small 
grid  of  V  =  1024£  sites.  The  six  plots  on  the  right  show  the  individual  particle  occupancies  of  the  microscopic 
state  at  time  step  intervals  of  128r.  The  three  local  states  per  site  of  the  lD3Px  model  are  plotted  in  the  vertical, 
respectively  labeled  bits  0,  1,  and  2.  A  black  line  marked  on  a  local  state  indicates  that  that  local  state  is  occupied 
by  a  single  particle.  The  profiles  on  the  left  show  the  mass  density  obtained  by  coarse-grain  averaging  over  blocks  of 
128  microscopic  cells  ( i.e .  total  mass  of  all  particles  in  a  block  divided  by  the  block  size).  On  the  upper  left  plot,  the 
red  sinusoid  is  the  initial  profile  at  t  =  Or  and  subsequent  profiles  at  t  =  12 8r  and  t  =  256r  are  overplotted.  On  the 
lower  left  plot,  the  profiles  at  t  =  384,  512,  and  640r  are  shown.  The  final  sinusoid  occurring  at  t  =  640r  is  rendered 
in  blue. 


Here,  the  angled  brackets  around  a  microscopic  quantity  denote  its  mesoscopic  expectation  value 
obtained  by  ensemble  averaging.  The  kinetic  transport  equations  are 

fo(x,t  +  r)  =  f0(x,t)  +  (n(x,t))  (3.9) 

f1>2(x±£,t  +  T)  =  f1>2{x,t)  - 

To  carry  out  a  classical  lattice-gas  simulation  at  the  mesoscopic  scale,  we  can  approximate  fimeso(x,  t)  = 
(f2(x,  i))  by  a  mean-field  collision  term,  denoted  fimf(x,f),  that  neglect  particle-particle  correlations 

<fi(x,t)>  ~  =  ft  h  (i  -  fo)  A)(i  -  /2).  (3.io) 

A  statement  of  detailed  balance  can  be  written  by  setting  the  mean-field  value  of  the  collision  term, 
Equation  (3.10),  to  zero  at  equilibrium 


(fi)  ~  nmf(/*eq)  =  0.  (3.11) 

Therefore,  the  probability  of  occupancies  satisfies  the  following  equation 

req  re q 

£ecL  _  _ J  1  J  2 _  /o  1  p\ 

Jo  ~  /r/2q  +  (i-/r)(i-/2eq)' 

This  equation,  along  with  equations  for  the  mass  and  momentum  densities 

Po  =  2/0-  +  /r  +  /2-  and  uXo  =  /r-/2eq,  (3.13) 

gives  us  a  nonlinear  system  of  three  equations  in  five  unknowns,  /^q,  /|q,  pa,  and  uxo.  Hence,  it 

is  possible  to  analytically  solve  for  the  occupation  probabilities,  /)5q,  and  fp,  in  terms  of  pQ  and 
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Figure  3.3:  A  bounding  polytope  for  the  lD3Px  lattice  gas  is  plotted  in  the  upper  left  corner  (here  the  mass  on 
the  x-axis  is  in  units  of  the  individual  particle  mass,  denoted  m).  The  other  plots  are  the  computed  equilibrium 
probability  occupancies.  The  units  for  momentum  on  all  the  plots  are  m^. 
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uxo.  The  resulting  functions  are  plotted  in  Figure  3.3.  The  shaded  region  region  in  the  upper  left  plot 
gives  the  allowable  values  for  mass  and  momentum  for  the  lD3Px  lattice-model.  The  equilibrium 
probability  of  occupancy  of  the  rest  particle,  /gq,  and  the  speed  ±1  moving  particle  states,  /qq  and 
fP,  can  be  analytically  determined  for  the  lD3Px  model.  The  results  are  plotted  for  a  range  of 
masses  from  m  =  0.5  up  to  m  =  3.5.  The  particle-hole  symmetry  about  half-filling  {in  =  2)  of  this 
lattice-gas  model  is  apparent  in  Figure  3.3  (for  example,  the  m  =  0.5  plot  is  identically  inverted  from 
the  m  =  3.5  plot).  fp  is  concave  downward  below  half-filling  and  concave  upward  above  half-filling. 
Furthermore,  below  half-filling,  fp  monotonically  increases  with  increasing  momentum  while  fp 
monotonically  decreases.  The  situation  is  reversed  above  half-filling,  as  expected. 

When  the  system  is  at  rest  at  equilibrium,  px  =  0,  then  fp  =  fp  =  d  and  the  probability  of 
occupancy  for  the  rest  particle  state  is 


feq  _ 

Jo  — 


1  -  2d  +  2d2' 


(3.14) 


Using  Equation  (3.10),  then  the  Jacobian  of  the  collision,  Ja & 


“  ~BfT 


/  — 1  +  2d  —  2d: 

1-2  d  +  2d2 
\  1  —  2d  +  2d2 

(1 —d)d 

(1  -d)d 

J  = 

l-2d+2d2 

(d—X)d 

l-2d+2d2 

(d-l)d 

l-2d+2d2 
( d-l)d 

l-2d+2d2 

{d-l)d 

l-2d+2d2 

l-2d+2d2 

The  eigenvectors  of  J  are 

|1> 

=  (2,1,1) 

|2> 

=  (0,1, 

-1) 

|3> 

-  2d  +  2d2)2 
d(d  —  1) 

,l,l)  ■ 

(3.15) 


(3.16) 

(3.17) 

(3.18) 


The  eigenvectors  |1)  and  |2),  corresponding  to  mass  and  momentum,  span  a  two-dimensional  hy¬ 
drodynamic  subspace.  The  remaining  eigenvector,  1 3) ,  is  a  kinetic  eigenvector,  which  in  this  case  is 
density  dependent.  The  eigenvalues  of  J  are 


Ai  =  0 

A2  =  0 

1  -  2d  +  6d2  -  8d3  +  4 d4 
3  “  -1  +  2 d-  2d2 


(3.19) 

(3.20) 

(3.21) 


Now  using  the  lattice  vectors,  eo  =  0,  ei  =  1,  and  e2  =  —  1,  and  the  expression  for  J  given  in 
Equation  (3.15),  setting  the  determinant  the  linearized  lattice-Boltzmann  equation  equal  to  zero 


^gi(faeaT+ajr) 


1)1- J 


=  0 


(3.22) 


allows  us  to  solve  to  the  dispersion  relations  for  the  system.  Taking  £  =  r  =  1,  we  get  the  following 
dispersion  relation 


(1-2 d+  2d2)e3u  -  2[d  -  3 d2  +  4 d3  -  2d4  +  (1  -  3 d+  3d2)  cos  k\e2uj+  (3.23) 

(1  -  2d)2[l  +  2 d{d  -  1)  cos  k\eu  +  4 d2{d  -  l)2  =  0.  (3.24) 


This  is  a  cubic  equation  in  eu ,  and  it  is  analytically  solvable.  The  only  hydrodynamic  mode  is  a 
damped  sound  wave,  as  expected.  Plots  of  the  dispersion  relation  for  various  background  densities 
are  given  in  Figure  3.4. 

Real  and  imaginary  parts  of  the  dispersion  relations  for  the  lD3Px  lattice-gas  model,  shown 
respectively  on  the  left  and  right  hand  sides  of  the  Figure  3.4.  The  real  part  of  the  dispersion 
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Figure  3.4:  The  real  and  imaginary  parts  of  the  dispersion  relations  for  the  lD3Px  lattice  gas  for  a  range  of 

background  densities  from  d  =  0.05  to  d  =  0.2. 
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relations  indicates  a  sound  mode  (Sft(u;)  — >  ±csk  as  k  — >  0).  The  imaginary  part  of  the  dispersion 
relation  for  the  hydrodynamic  mode  is  parabolic  for  small  wavenumber,  indicating  viscous  damping 
of  the  sound  mode  (S(w)  — >  Tk2  as  k  — >  0).  Remarkably,  the  sound  damping  constant,  T,  approaches 
zero  as  the  background  mass  density  approaches  zero.  That  is,  low-mass  density  waves  can  oscillate 
without  viscous  damping.  This  behavior  is  observed  in  numerical  simulations  of  the  lD3Px  model. 


0  1  2  3  4  5  6 


Wave  Number  (k) 


Wave  Number  (k) 


Figure  3.5:  Real  part  of  the  dispersion  relation  for  the  mesoscopic  lD3Px  lattice  gas  in  the  long  wavelength  limit 
and  mean-field  limit  at  a  reduced  background  density  of  d  =  0.214286. 


The  real  part  of  the  dispersion  relation  for  the  sound  mode  for  the  lD3Px  lattice-gas  model  set 
with  a  background  density  of  d  =  with  V  =  7,  is  shown  in  Figure  3.5.  The  real  part  of  the 
dispersion  relation  indicates  a  sound  mode  (3?(w)  — >  ±csk  as  k  — >  0  where  cs  =  0.74|).  The  data 
points,  plotted  as  black  circles,  are  solutions  to  the  linearized  Boltzmann  equation  in  the  mean- field 
limit.  The  solid  red  curves,  with  slope  of  ±cs,  are  numerical  linear  fits  to  the  data.  The  imaginary 
part  of  the  dispersion  relation  for  the  sound  mode  for  the  lD3Px  lattice-gas  model  is  shown  in 
Figure  3.6.  The  imaginary  part  of  the  dispersion  relation  indicates  sound  damping  (S(w)  — >  iTk 2  as 
k  — >  0  where  T  =  0.08 ^r.  The  solid  red  parabola  is  a  numerical  fit  to  the  data  in  the  region  of  small 
k  <  1.  The  calculations  shown  in  Figures  3.5  and  3.6  were  done  with  a  mass  density  filling  fraction 
of  d0  =  ^  =  0.214,  where  a  small  system  size  of  V  =  7  is  used.  In  this  case,  k  =  -y-  =  0.898. 


3.3  Comparing  Analytical  and  Numerical  Predictions 

A  time  history  of  the  mass  density  wave  for  a  small  system  with  V  =  7  sites  is  shown  in  Figure  3.7. 
The  exponential  envelope  over  plotted  in  red  analytically  determined  by  an  analysis  of  the  linearized 
lattice  Boltzmann  equation  in  the  mean-field  limit  (see  Figure  3.6).  The  predicted  sound  damping 
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Figure  3.6:  The  imaginary  part  of  the  dispersion  relation  for  the  mesoscopic  lD3Px  lattice  gas  in  the  long 

wavelength  limit  and  mean-field  limit  at  a  reduced  background  density  of  d  =  0.214286. 


Time 


Figure  3.7:  Damping  of  a  mass  density  wave  for  a  system  with  V  —  7  sites  in  the  classical  lD3Px  model  simulated 
using  a  mesoscopic  Boltzmann  equation  with  the  collision  term  expressed  in  the  mean-field  approximation.  The 
background  density  is  d0  =  =  0.214.  The  ordinate  is  the  absolute  value  of  the  amplitude  of  mass-density  wave 

divided  by  the  peak  amplitude  of  the  initial  perturbation. 
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p2 

constant,  T  =  0.08^  is  in  excellent  agreement  with  the 


simulation  data. 


Time 


Time 


Figure  3.8:  Damping  of  a  mass  density  wave  at  multiple  scales  in  the  classical  lD3Px  model  simulated  using  a 
mesoscopic  Boltzmann  equation  with  the  collision  term  expressed  in  the  mean-field  approximation.  The  abscissa  is 
plotted  in  units  of  r  and  the  ordinate  is  plotted  as  the  absolute  value  of  the  amplitude  of  mass-density  wave  divided 
by  the  peak  amplitude  of  the  initial  perturbation. 

Damping  of  mass  density  waves  in  the  classical  lD3Px  lattice  gas  for  different  system  sizes, 
V  =  2,3,4, 5,8  and  16  is  shown  in  Figure  3.8.  The  simulation  is  initialized  with  a  sinusoidal 
perturbation  of  Sp  =  0.04  from  a  uniform  background  mass  density  at  half-filling,  p  =  2.  So  the 
fractional  mass  density  variation  is  initially  one  part  in  fifty.  The  wavelength  equals  the  system 
size.  Plotted  is  the  time  history  of  the  resulting  standing  wave.  The  peak  amplitude  decays  in  time 
because  of  viscous  damping.  The  circles  are  data  taken  from  classical  mesoscopic  simulations,  using 
a  mean-field  collision  operator.  The  solid  blue  curve  is  a  nonlinear  numerical  fit  to  the  data  using 
an  exponentially  damped  sinusoid  e~  ?  cos  cot  with  free  parameters  for  the  angular  frequency,  to,  and 
the  damping  constant,  r.  The  dotted  blue  curve  is  the  envelope  of  the  damped  curve.  For  system 
sizes  as  small  as  V  =  3,  the  oscillations  are  clearly  apparent.  Even  in  the  smallest  possible  system 
only  with  two  sites,  V  =  2,  oscillations  are  discernable,  and  the  data  agrees  with  an  exponentially 
damped  sinusoid.  This  is  an  example  of  “fluid- like”  behavior  occurring  in  systems  far  below  from 
the  continuum  limit. 

Plotted  in  Figure  3.9  is  damping  time  constants  of  mass  density  waves  in  the  classical  lD3Px 
lattice  gas  for  different  system  sizes  from  V  =  2  up  to  V  =  256.  The  log-log  plot  shows  the  power-law 
behavior,  known  as  diffusive  ordering ,  typical  of  lattice-gas  system  in  the  viscous  regime.  The  power 
law  in  this  case  is  T  =  0.44V2,  which  is  parabolic.  Each  circle  is  determined  from  a  mesoscopic  scale 
simulation  that  initialized  with  a  sinusoidal  perturbation  of  Sp  =  0.04  from  a  uniform  background 
mass  density  at  half-filling,  p  =  2.  The  damping  constant,  T  =  y,  is  determined  from  the  envelope 
of  the  resulting  standing  wave  e~i  cos  cot  (see  Figure  3.8).  The  mean-field  estimates  of  the  damping 
time  constant  are  the  circles.  The  solid  blue  line  is  a  linear  best  fit  to  these  estimates.  The  estimated 
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Figure  3.9:  Diffusive  ordering  in  the  classical  lD3Px  model  computing  at  the  mesoscopic  scale  using  the  mean-field 
approximation. 


damping  constant  deviates  only  slightly  from  power-law  behavior  at  the  smallest  system  sizes.  The 
inset  plot  is  a  linear  plot  of  the  data  for  V  <  16  and  the  solid-blue  parabola  is  the  same  diffusive- 
ordering  power-law  in  the  larger  log-log  plot. 


Chapter  4 

A  Two-Dimensional  Model  with 
Conserved  Mass  and  Momentum 


4.1  Microscopic  Rules 


In  the  two-dimensional  lattice  gas  on  a  triangular  lattice  ( B  =  6)  there  are  9  equivalence  classes: 

•  a  three  member  equivalence  class  for  the  2-body  head-on  collisions  p  =  0  ; 

•  a  two  member  equivalence  class  for  the  symmetric  3-body  collisions  with  p  =  0; 

•  six  two  member  equivalence  classes  for  2-body  head-on  collisions  with  a  spectator  particle  (or 
3-body  asymmetric  collisions)  with  momentum  p  =  mcea ;  and,  finally, 

•  a  three  member  equivalence  class  for  the  4-body  collision  with  p  =  0,  which  is  the  particle-hole 
counterpart  of  the  2-body  equivalence  class. 


The  equivalence  classes  are  illustrated  in  Figure  4.1.  The  first  two  are  the  equivalence  classes  used  in 
the  FHP  model  [2].  There  are  three  conserved  quantities  for  this  two-dimensional  system:  the  mass, 
and  two  components  of  the  momentum.1  Energy  is  also  conserved,  but  this  is  degenerate  with  the 
mass  since  the  particles  move  at  unit  speed.  The  FHP  model  does  not  use  the  other  7  equivalence 
classes,  even  though  adding  them  in  takes  no  extra  computational  or  analytical  work. 

Let  the  incoming  configuration  at  a  spacetime  coordinate  ( x ,  t)  be  denoted  by  the  set  s(x,  t)  = 
{ni(x,  f),  ri2(x,  t), . . . ,  ns{x,  t)}.  After  the  collision  step,  the  outgoing  configuration  is  denoted 
s'(x,t)  =  {n\  (x.  t) ,  rig (x,  t) , . . . ,  n’B  {x.  t) } .  Taking  B  =  6  for  example,  a  3-body  collision  can  have 
an  incoming  configuration 

{10  10  10} 


and  an  outgoing  configuration 


{010101}, 


or  vice  versa.  In  the  case  of  a  2-body  collision,  if  the  incoming  configuration  is 


{1  0  0  1  0  0}, 


then  the  outgoing  configuration  could  be  either 


{010010}  or  {0  0  1  0  0  1}. 

1  If  we  use  only  the  2-body  collisions  in  Figure  4.1,  then  there  is  would  be  additional  conserved  quantity  (the 
difference  in  the  particle  number  along  each  of  the  three  lattice  directions  give  three  conserved  momenta  instead  of 
two).  So  there  would  be  a  spurious  invariant.  Consequently,  the  symmetric  3-body  collisions  in  Figure  4.1  are  needed 
in  the  FHP-model  to  remove  this  spurious  invariant  [2]. 
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Figure  4.1:  Successive  snap  shots  of  the  time  evolution  of  a  small  classical  lattice-gas  system  on  a  two-dimensional 
8x8  triangular  grid.  The  particles  are  rendered  in  red  and  the  underlying  hexagonal  cells  of  the  Bravais  lattice  are 
rendered  in  blue.  Each  hexagonal  cell  comprises  6  local  states  and  hence  each  cell  can  hold  up  to  6  particles.  There 
are  nine  collision  sets,  called  equivalence  classes.  The  first  three  head-on  m  =  2  collisions  and  two  symmetric  m  —  3 
collisions,  shown  in  the  red  box,  constitute  the  FHP  lattice-gas  model.  (Viscous  dissipation  can  be  reduced  if  the 
m  =  3  asymmetric  collisions  and  the  m  =  4  collisions  are  used,  but  they  are  not  necessary.)  Initially,  at  t  —  0,  most 
of  the  particles  reside  on  the  left  half  of  the  lattice.  (The  sinusoidal  profile  across  the  lattice  is  highly  exaggerated 
for  display  purposes,  so  the  system  is  very  far  from  local  equilibrium.)  After  a  few  time  steps,  at  t  =  5,  the  situation 
is  reversed  as  most  of  the  particles  have  moved  to  local  states  on  the  right  half  of  the  lattice.  Each  configuration  of 
particles  rendered  in  this  figure  corresponds  to  a  unique  state  of  the  lattice-gas  system.  Even  though  this  is  a  small 
system  with  only  64  lattice  sites,  the  total  number  of  unique  states  available  is  still  a  rather  large  number,  2384.  The 
CAM-8,  our  fastest  lattice-gas  computer,  which  can  update  25  million  sites  per  second,  would  have  to  run  for  over 
10lo°  years  to  enumerate  every  possible  state  of  this  small  8x8  example. 
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The  collision  part  of  the  classical  dynamics  permutes  the  particles  locally  because  the  total  number 
of  particles  and  the  total  momentum  must  be  conserved.  The  collision  equation,  which  is  applied 
homogeneously  across  the  lattice,  is  the  following 


s'  =  U(s), 


lookup  table  defined  by 


where  U  is  a  64  element 

f7({100100}) 
[/({010010}) 
f7({001001}) 
t/  ({101010}) 
[/({010010}) 


=  a{010010}  +  (1  -  a){001001} 
=  a{001001}  +  (l-a){100100} 
=  a{100100}  +  (1  -  a){010010} 
=  {010101} 

=  {101010}, 


where  the  a  €  {0, 1}  is  chosen  uniformly  randomly.  If  a  configuration  s  is  not  a  member  of  the 
either  the  2-body  or  the  3-body  FHP  equivalence  classes,  then  U(s)  =  s.  In  other  words,  if  the 
incoming  state  is  not  a  member  of  an  equivalence  class,  which  has  two  or  more  members,  then  the 
outgoing  state  is  set  equal  to  the  incoming  state.  Some  fluid  simulation  examples  of  the  FHP  model 
are  shown  in  Figures  4.2  and  4.3. 


Figure  4.2:  Classical  lattice  gas  fluid  simulation  with  a  lattice  size  of  4096  X  2048  done  on  the  CAM-8  using  about 
10  million  particles.  Momentum  and  vorticity  fields  display  a  Von  Karman  Street.  Coarse-gain  averaging  was  done 
over  50  time  steps,  using  a  64  X  64  spatial  block  size  for  the  momentum  field,  and  a  16  X  16  spatial  block  size  for  the 
vorticity  field.  Red  indicates  clockwise  vorticity  and  blue  counter  clockwise  vorticity.  The  system  was  initialized  with 
a  reduced  mass  density  of  d  =  1/7.  The  system  was  at  rest  at  t  =  0,  and  was  accelerated  to  the  right  by  external 
forcing  resulting  in  a  terminal  flow  velocity  of  v  =  0.3 c  at  t  =  20,  000,  shown  in  this  figure.  The  diameter  of  the 
cylindrical  obstacle  is  256 i  lattice  units.  The  critical  Reynolds  number  for  vortex  shedding,  Re  =  42,  was  achieve  at 
approximately  t  =  10,  000  time  steps  into  the  simulation.  The  maximum  Reynolds  number  achieved  in  simulation  is 
about  Re  ~  250. 


4.2  Model  Analysis 

The  FHP  collision  term  is 

=  H^=2.p=o  +  (4.1) 

The  two-body  term  is 

^  —  ^^a+l^a+4  (1  ^a)(l  ^a+2)(l  ^a+3)(l  ^a+ 5) 
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Figure  4.3:  Classical  lattice-gas  fluid  simulation  with  about  10  million  particles  computed  on  the  CAM-8.  Vorticity 
and  momentum  fields  of  the  two-dimensional  shear  instability  are  shown.  The  simulation  was  done  using  a  lattice  size 
of  4096  X  2048  with  toroidal  boundary  conditions,  with  spacetime  averaging  over  128x128  blocks  for  50  time  steps. 
FHP  collisions  with  spectators  and  a  rest  particle  were  implemented.  Macroscopic  scale  data  is  presented  at  0,  10000, 
and  30000  time  steps,  which  is  at  a  Reynolds  number  of  approximately  Re  ~  1000. 
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+  (1  -  a)na+2na+5(l  -  na)(  1  -  na+ i)(l  -  na+3)(l  -  na+4) 

-  «a«a+3(l  -  ^a+i)(l  -  na+2){l  -  na+4)(l  -  na+5),  (4.2) 


where  a  £  {0,1}  is  chosen  uniformly  randomly  and  and  the  symmetric  3-body  term  is 

0.a  ^  —  ria+ina+3?rQ+5(l  na)(l  ?ra4-2)(l  ^a-t-4) 

-  nQ7iQ+2na+4(l  -  na+i)(l  -  na+3)(l  -  na+5).  (4.3) 


The  indices  are  taken  modulo  B.  Writing  fiPHP(?i*)  with  and  asterisk  subscript  on  n*  denotes  that 
the  collision  term  for  the  ath  local  state  depends  on  all  the  on-site  local  states. 

The  Jacobian  of  the  mesoscopic  mean- field  collision  term  for  the  FHP  lattice  gas  is  a  circulant 
matrix2 


JFHP 


/  —d(l—d)2 

±d(l+d){l-d)2 
ld(l-3d)(l-d)2 
-d(l-2d)(l-d)2 
id(l-3d)(l-d)2 

V  id(l+d)(l-d)2 


±d(l+d)(l-d)2 

-d(l-d)2 

id(l+d)(l-d)2 

id(l-3d)(l-d)2 

-d(l-2d)(l-d)2 

id(l-3d)(l-d)2 


id(l-3d)(l-d)2 
id(l+d)(l-d)2 
-d(l-d)2 
id(l+d)(l-d)2 
id(l-3d)(l-d)2 
—  d(l  — 2d)(l  — d)2 


—  d(l  — 2d)(l  — d)2 
id(l-3d)(l-d)2 
id(l+d)(l-d)2 
-d(l-d)2 
id(l+d)(l-d)2 
}d(l-3d)(l-d)2 


id(l-3d)(l-d)2 
—  d(l  — 2d)(l  — d)2 
id(l-3d)(l-d)2 
id(l+d)(l-d)2 
-d(l-d)2 
id(l+d)(l-d)2 


id(l+d)(l-d)2  \ 
id(l-3d)(l-d)2 
—  d(l  — 2d)(l  — d)2 
id(l-3d)(l-d)2 
}d(l+d)(l-d)2 
-d(l-d)2  / 

(4.4) 
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Figure  4.4:  Dispersion  relations  for  the  FHP  lattice  gas  for  k  along  the  rr-axis  direction,  u  =  uj(p,kx)  at  a 

background  density  of  p  =  2.5.  The  real  part  of  uj(kx ),  plotted  on  the  left,  indicates  propagating  modes,  and  the 
imaginary  part  of  tu(kx),  plotted  on  the  right,  indicates  two  damped  soft  excitations  with  parabolic  dispersion  relations 
at  kx  ~  0  and  spurious  modes  at  finite  kx.  The  periodicity  of  the  reciprocal  lattice  is  observed  in  this  numerical 
solution,  for  example  inversion  symmetry  u(—k)  =  uj(k). 


For  the  FHP  lattice  gas  there  are  three  soft  modes,  one  corresponding  to  conservation  mass  and 
two  corresponding  to  the  two  conserved  components  of  momentum.  The  kinetic  modes  which  are 
nonvanishing  at  first  order  cause  damping  in  the  lattice  gas  and  are  attributed  to  a  positive  imaginary 
part  of  a;,  (Im(a;)  >  0).  There  are  three  hard  modes  corresponding  the  transport  coefficients  of  the 
bulk,  shear,  and  cubic  viscosities,  and  in  the  long  wavelength  limit,  the  shear  and  cubic  viscosities 
coincide  for  an  isotropic  fluid. 

2  It  is  sometimes  possible  to  go  beyond  this  Boltzmann  approximation  to  determine  a  more  accurate  form  of  J 
by  accounting  for  the  particle-particle  correlations,  which  arise  from  on-site  collisions  at  the  microscopic  scale.  This 
is  called  renormalization ,  and  it  which  give  rise  to  corrections  the  Boltzmann  estimates  of  the  transport  coefficients. 
Renormalization  theory  for  lattice  gases  has  been  worked  out  by  Ernst  et  al.  [29,  31]  and  by  Boghosian  [88]  who  has 
developed  diagrammatic  expansion  methods  for  determining  a  renormalized  J.  Various  approximations  to  sum  terms 
involving  higher  and  higher  order  particle-particle  correlations,  including  the  ring  collisions  for  example,  in  certain 
cases  are  necessary  to  correctly  estimate  the  value  of  the  transport  coefficients  of  a  lattice-gas  system.  However,  the 
renormalization  of  J  is  not  necessary  for  single-speed  lattice  gases  where  the  Boltzmann  estimates  are  quite  accurate. 
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Figures  4.5  shows  a  smaller  region  of  the  k- space  within  the  first  Brillouin  zone  of  the  reciprocal 
lattice  where  the  long  wavelength  behavior  of  the  system  is  apparent.  A  sound  mode  is  observed  in 
the  real  part  of  u>.  For  k  along  x  and  y- axis,  the  dispersion  relation  in  the  long  wavelength  limit 
predicts  a  sound  speed  of  cs  =  ^,  which  agrees  with  the  numerical  measurements.  Since  cs  is  the 
same  when  the  wave  vector  is  aligned  with  either  the  x  or  y- axes,  this  indicates  that  the  sound  mode 
is  isotropic,  see  the  L.H.S.  of  Figures  4.5  and  4.6.  The  imaginary  part  of  w  is  parabolic  at  k  — >  0. 
This  is  characteristic  of  a  viscous  hydrodynamics  regime  with  diffusive  ordering,  see  the  R.H.S.  of 
Figures  4.5  and  4.6.  In  the  viscous  hydrodynamic  regime,  the  dispersion  relations  for  sound  mode 
in  the  FHP  lattice  gas  is  found  to  be 


uj(k)  =  ±csk  +  iT{p)k2 .  (4-5) 

Higher  resolution  plots  of  the  dispersion  curves,  for  a  large  range  of  background  mass  densities,  are 
shown  in  Figure  4.13.  We  numerically  fit  the  shallowest  mode  with  a  parabola.  This  gives  us  a  way 
to  determine  the  damping  constant  of  compressional  mass  density  waves.  This  result  is  given  in 
§4.3.3 


0.2  0.4  0.6  0.8  1 


Figure  4.5:  Dispersion  relations  for  the  FHP  lattice  gas  for  k  along  the  £-axis  direction,  w  =  uj(p,kx)  at  a 

background  density  of  p  =  2.5  and  for  kx  <  1.  In  the  long  wavelength  limit  (k  —t  0) ,  the  real  part  of  uj(kx),  plotted  on 
the  left,  predicts  a  propagating  excitation  which  is  the  sound  mode  with  a  linear  dispersion  relation  SR-fiuf/Ca,)}  =  ±esfcx 
with  sound  speed  cs  =  -y=.  The  imaginary  part  of  uj(kx),  plotted  on  the  right,  indicates  two  damped  soft  excitations 
with  parabolic  dispersion  relations  at  kx  — >  0.  The  lowest  parabola  indicates  sound  damping  with  the  dispersion 
relation,  ^s{uj(kx)}  =  T/c^,  which  is  plotted  in  more  detail  and  for  a  range  of  densities  in  Figure  4.11. 


Here  we  illustrate  the  identity  (2.24) 

J~X\eiej)  =  T~leiei)  (4-6) 

I^Tj 

using  the  well-known  two-dimensional  FHP  lattice-gas  model  [2].  Using  the  x  and  y  components  of 
ea  =  (cos  ^,sin  separately,  we  define  the  following  two  kets 

\ex)  = 

\ev)  = 

The  direct  product  of  |ej)  and  | ej )  for  (i.j)  G  {x,y)  component  by  component  gives 

/q 

\exey)  =  ^(1,  -1,°, -1,1,0) 


*(1,-1,  -2, -1,1,  2) 
^(1,1, 0,-1, -1,0). 


(4.7) 
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Figure  4.6:  Dispersion  relations  for  the  FHP  lattice  gas  for  k  along  the  y-axis  direction,  cj  =  uj(p,ky)  at  a 

background  density  of  p  =  2.5.  In  the  long  wavelength  limit  ( k  — »•  0),  the  real  part  of  uj(ky),  plotted  on  the  left,  predicts 
a  propagating  excitation  which  is  the  sound  mode  with  a  linear  dispersion  relation  =  zb csky  =  cs27r-^F- 

with  sound  speed  cs  =  ~^=.  The  factor  of  sin(^)  =  comes  about  because  the  triangular  lattice’s  cell  size  is 
shorter  than  it  is  wide  by  this  amount.  The  imaginary  part  of  to(kx),  plotted  on  the  right,  indicates  two  damped  soft 
excitations  with  parabolic  dispersion  relations  at  kx  — >  0.  The  lowest  parabola  indicates  sound  damping  with  the 
dispersion  relation,  S s{u(ky )}  =  Fky.  At  ky  0,  this  is  identical  to  the  dispersion  relation  with  wave  vector  along 
the  x-axis  indicating  an  isotropic  sound  mode. 


\e-xex)  =  ^(1,1,4, 1,1,4) 


\evev)  ~ 


73 


(1, 1,0, 1,1,0). 


(4.8) 

(4.9) 


Following  Wolfram  [1],  we  take  the  basic  approach  that  the  eigenkets  of  the  linearized  collision 
operator,  which  is  the  circulant  matrix  for  the  FHP  model,  have  the  form 


y)  = 


t  e  T„  \ 

— e  1  3 
(-1)“ 

— e'T 


(4.10) 


for  a=l,2,...,6.  We  are  free  to  define  a  new  set  of  eigenkets  by  taking  linear  combinations  of  |a) 
so  that  we  can  express  the  eigenkets  of  JFHP  in  terms  of  the  lattice  kets  \ex),  |ey),  le^e^),  \eyey), 
and  \exey).  This  is  done  as  follows 

|1')  =  |6)  =  (l,l,l,l,l,l)  =  |exex)  +  |e^) 

|2'>  =  ^(|l)-|5))  =  (l,l,0,-l,-l,0)  =  -T|ey) 

|3'>  =  |1)  +  |5)  =  (1,  —1,  —2,  —1, 1, 2)  =  2|ex) 

|4')  =  |2)  —  |4)  =  (1, 1,  —2, 1, 1,  —2)  =  2(\eyey)  —  |exex» 

|5')  =  ^(|2)  +  |4))  =  (l,-l,0,l,-l,0)  =  -T|exey) 

|6')  =  |3)  =  (—1, 1,  —1, 1,  —1, 1). 


All  the  ( a'\  must  be  defined  with  the  appropriate  magnitude  to  satisfy  the  normality  condition, 
( a'\a ')  =  1  for  all  a'.  The  mean-field  collision  operator  for  the  FHP  lattice  gas  is  simply  obtained 
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by  replacing  the  na’s  in  Equation  (4.1)  with  /a’s.  Using  this  in  the  eigenvalue  formula  (2.18) 

B 

Ka=J2  Jme2~/B,  (4.11) 

m—  1 

it  follows  that  the  eigenvalues  of  JFHP  are 


(!' 

|JPHP 

|1') 

=  0 

(2' 

|JPHP 

|2') 

=  0 

(3' 

|JPHP 

13') 

=  0 

(4' 

|JPHP 

|4') 

=  — 3d(l  - 

- d )3 

(5' 

|JPHP 

15') 

=  — 3d(l  - 

- d )3 

(6' 

|JPHP 

|6') 

=  6  c?2  (1  - 

df. 

Note  that  (l'|  JFHP|1/)  corresponds  to  mass  conservation,  (2'|  JFHP|2')  to  y-momentum  conservation, 
and  (3'|  JFHP|3')  to  x-momentum  conservation.  The  degenerate  viscous  eigenvalue  kv  =  — 3d(l  —  d )3 
is  immediately  identified. 

Using  the  eigenket  of  the  kinetic  subspace,  we  can  compute  the  generalized  inverse  (JPHP)_1  as 
follows 


(vT 

0 

°\ 

/(4'l\ 

( JFHP)_1  =  ( |4')  5')  6') ) 

0 

1 

0 

(5'| 

S - V - ' 

(6x3  matrix) 

l 0 

0 

i) 

V<6'iy 

(3x6  matrix) 


(4.12) 


We  have 


(JPHP)-1 


/  1+3  d 

-(1+d) 

1-3  d 

-1+5  d 

1-3  d 

—  (1+d) 

36  (  —  1  +  d)3  d2 

36  (  —  1  +  d)3  d2 

36  (  —  1+d)3  d2 

36  (  —  1+d)3  d2 

36  ( — 1+d)3  d2 

36  ( — 1+d)3  d2 

—  (1+d) 
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-(1+d) 
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1-3  d 
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1-3  d 

-(1+d) 
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—  (1+d) 

1-3  d 

36  (-1+d)3  d2 

36  (  —  1  — |—  tZ)  ^  d2 

36  (  —  1+d)3  d2 

36  (  —  1+d)3  d2 

36  ( — 1+d)3  d2 

36  ( — 1+d)3  d2 
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36  (  —  1  +  d)3  d2 
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36  (  —  1+d)3  d2 

36  ( — 1+d)3  d2 

36  (—1+d)3  d2 

—  (1+d) 

1-3  d 

-1+5  d 

1-3  d 

—  (1+d) 

1+3  d 

\  36  (-1+d)3  d2 

36  (  —  1  — |—  cZ)  ^  d2 

36  (  —  1+d)3  d2 

36  (  —  1+d)3  d2 

36  ( — 1+d)3  d2 

36  ( — 1+d)3  d2 

(4.13) 


Note  that  (JPHP)  1 JPHP  ^  l,  since  (</PHP)  1  is  the  generalized  inverse  of  JPHP 


/  2 


(JPHP)-1JFHP  = 


V-l 


h 


o  -k 


o 
1 
6 
0 

_  1 

7/ 

O  / 


(4.14) 


However,  if  we  compute  the  matrix  of  elements,  we  see  that  J  1  is  indeed  the  inverse  of  J  in  the 
kinetic  space 

/0  0  0  0  0  0\ 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  1  0  0 

0  0  0  0  1  0 

\o  0  0  0  0  lJ 


(a'  |  ( JFHP)  ” 1 JPHP  |  a')  = 


(4.15) 
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Finally,  we  see  that  \exev)  is  an  eigenvector  of  (JPHP)  1 ,  with  eigenvalue  ^  =  3d(1id)3 


{rn-'M  = 


*)■ 


v/  ~  3d(i~dr]  x  y/ 

The  Henon’s  shear  viscosity  formula  (2.50)  for  a  classical  lattice  gas  is 


V  = 


pi'i  1  (  1  1 


r  D  +  2  \  k,  2 


Since  nv  =  — 3d(l  —  d)3  for  the  FHP  lattice  gas,  Equation  (2.50)  is 


FHP 

V  =  p- 


t  \  12d(l  —  d)3  8 J 


(4.16) 


(4.17) 


(4.18) 


4.3  Comparing  Analytical  and  Numerical  Predictions 

4.3.1  Single-Particle  Probability  of  Occupancy  in  the  Subsonic  Limit 
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Figure  4.7:  Theory  versus  simulation  comparison  of  the  velocity  dependence  of  the  single-particle  distribution 
function  in  the  non-Galilean  parameterization:  fa  =  d  +  dDea  ■  v  -)-  gdD(D/2  +  1  )Qa  :  vv.  FHP  simulation  data 
is  overplotted  on  this  predicted  mesoscopic  probability  of  occupancy.  Plots  (a)  and  (b)  are  for  reduced  background 
densities  of  d  =  .20  and  d  =  0.25,  respectively.  A  velocity  shift  is  imparted  along  the  rr-axis;  that  is,  along  the  f\ 
direction  indicated  in  the  figure.  Data  was  collected  from  a  128  x  128  classical  FHP  simulation  (crosses)  and  was 
coarse-grained  averaged  over  1600  time  steps  from  time  step  t  =  400  to  t  =  2000. 


The  general  form  of  the  single-particle  occupancy  probability,  appropriate  for  single  speed  lattice 
gases,  is  a  Fermi-Dirac  function  whose  argument  is  the  sum  of  scalar  collision  invariants,  ap  4-  (3ea  • 
P  +  l E  [3] 


feq  = 

J  a 


l 


X  _|_  g«P+/3ea-p+7f5  ' 


(4.19) 


Fundamentally,  this  arises  because  the  individual  classical  bits  representing  particles  satisfy  a  Pauli 
exclusion  principle.  By  Taylor  expanding  (4.19)  about  v  =  0  to  fourth  order  in  the  velocity  and 
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equating  the  zeroth,  first,  and  second  moments  of  fa  respectively,  the  parameters  a  and  j3  are 
determined  [3] .  The  result  of  the  low  Mach  number  expansion  of  fqq  is 


n 

~B 


nD  nD(D  +  2) 

^  jj  t'aiVi  T  9  2c2B  ^'ai^aj'ViVj  Q 


n(D  +  2) 
2  c2B 


v2  +  0{M3), 


where  the  density  dependent  Galilean  prefactor  is 


(4.20) 


9(d) 


D  1-2  d 

D  +  2  1  -  d  ' 


(4.21) 


The  predicted  functional  form  (4.20)  is  checked  in  the  subsonic  limit  against  numerical  data  taken 
from  a  simulation  of  a  classical  FHP  lattice  gas.  The  analytical  and  numerical  results  are  in  excellent 
agreement  and  are  shown  in  Figure  4.7. 

Using  p  =  mn  for  the  density  and  cs  =  for  the  sound  speed,  the  moments  of  lattice  gas 
distribution  are 


»E(/;)‘"  =  P  (4.22) 

a 

mc^2eai(f^ydeal  =  pvt  (4.23) 

a 

mc2J2eaieaj{faqYdeal  =  pc2s{l  -  g^)5lj  +  gpviVj.  (4.24) 

a 

Note  that  for  a  uniform  filling  of  states,  fa=d  for  all  directions  and  so  v  =  0,  then3 

p(x,  t)  =  mBd  (4-25) 

Vi(x,t )  =  0  (4.26) 

II  ij{x,t)  =  pc25ij.  (4.27) 


This  is  expected  since  when  ty  =  0,  the  momentum  flux  density  tensor  is  diagonal  and  scales  with 
the  pressure  (II ^  =  P8ij)  and  for  an  ideal  gas  P  =  pc2. 

4.3.2  Measuring  the  Shear  Viscosity  Using  the  Decay  of  a  Sinusoidal  Shear 
Wave 

Given  a  sinusoidal  perturbation  of  wavelength  A  of  a  fluid  one  can  straightforwardly  measure 
the  time  for  relaxation  to  an  equilibrium  state  where  the  fluid  is  at  rest.  This  method  was  used  to 
measure  the  shear  viscosity  of  the  FCHC  lattice  gas  by  Adler  et  al.  [22].  The  relevant  part  of  the 
Navier-Stokes  equation  is  the  time  dependent  term  and  the  momentum  diffusion  term 

( dt  -  vd2)px  =  0.  (4.28) 

This  has  the  solution 

Px  =  p0sin(ky)e~k2vt.  (4.29) 

Therefore,  the  decay  rate,  k2v ,  can  be  measured  to  determine  v  since  k  =  ^  is  known.  This  method 
is  easier  to  implement  on  the  CAM-8  than  the  square-wave  forcing  method  used  by  Kadanoff,  McNa¬ 
mara,  and  Zanetti  [20],  since  no  forcing  bits  or  rules  are  required  and  it  is  easy  to  generate  an  initial 
random  fluid  pattern  with  a  sinusoidal  perturbation.  Very  good  agreement  is  found  between  the 
mean-field  prediction  of  the  kinematic  shear  viscosity  and  the  numerical  data  shown  in  Figure  4.10 
taken  on  the  CAM-8  for  the  FHP  lattice  gas. 

3  We  used  the  following  property  of  the  displacement  vectors,  ^  ea  =  0  and  )T]  £aiSaj  =  that  was 

originally  derived  by  Wolfram  [1] 
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Figure  4.8:  A  depiction  of  an  intial  sinusoidal  shear  wave  perturbation  that  decays  over  time  to  a  final  state  with 
a  smaller  amplitude. 


Time  Step/ 100 


Figure  4.9:  Plot  of  the  amplitude  of  a  shear  wave  versus  time.  The  amplitude  of  shear  wave  of  the  momentum 
density  field  decays  exponentially  as  is  shown  here  on  a  log-linear  plot.  Linear  regression  is  used  to  fit  the  data  set 
and  the  predicted  the  decay  exponent,  which  is  proportional  to  the  kinematic  shear  viscosity  of  the  lattice-gas  fluid. 
Data  is  taken  from  a  512  x  512  classical  FHP  simulation  at  a  background  density  of  d0  =  0.15. 
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Figure  4.10:  Kinematic  viscosity  versus  density  obtained  by  measuring  the  rate  of  exponential  damping  of  a 

sinusoidal  velocity  perturbation.  The  theoretical  mean-field  prediction  and  numerical  data  are  plotted  for  an  FHP 
lattice  gas  with  2  and  3-body  collisions  on  a  two  dimensional  triangular  lattice.  Simulation  runs  were  done  on  the 
CAM-8  on  a  512  x  512  periodic  space. 
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4.3.3  Measuring  the  Bulk  Viscosity  Using  the  Decay  of  a  Sinusoidal  Compres- 
sional  Wave 

In  the  viscous  hydrodynamic  regime,  the  dispersion  relations  for  sound  mode  in  the  FHP  lattice 
gas  has  the  form 

u>(k)  =  ±csfc  +  iT(p)k2. 

We  can  determine  the  damping  constant  by  solving  for  the  roots  of  the  secular  determinant  of  the 
linearized  Boltzmann  equation  in  fc-space.  Our  predicted  value  of  the  damping  constant  are  gotten 
by  a  numerical  fitting  procedure  shown  in  Figure  4.11.  In  this  case,  the  analysis  was  carried  out 
with  the  wave  vector  of  the  density  perturbation  aligned  with  the  z-axis  of  the  Bravais  lattice.  The 
procedure  was  repeated  for  a  density  perturbation  along  the  y-axis,  and  the  same  result  was  obtained. 
This  indicates  that  the  sound  mode  of  the  FHP  lattice  gas  is  isotropic  in  the  long  wavelength  limit 
where  diffusive  ordering  holds. 
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Figure  4.11:  Imaginary  part  of  the  dispersion  relation,  uj  =  uj(p,kx)},  for  the  FHP  lattice  gas.  The  wave  vector  is 
directed  along  the  £-axis,  k  =  kxx ,  and  the  results  are  plotted  for  a  range  of  background  densities  from  p0  =  0.6  to 
p  =  3.  The  lowest  curves  are  numerically  fit  with  a  parabola,  which  are  plotted  in  red.  Since  the  dispersion  relation 
is  ^s{cu(ky )}  =  T k%.  in  the  long  wavelength  limit,  we  can  predict  the  damping  constant  as  a  function  of  mass  density, 
T  =  T(p).  This  analytical  prediction  of  T(p)  agrees  well  with  data  from  numerical  measurements  taken  from  the 
lattice-gas  simulation.  See  Figure  4.13. 


46  CHAPTER  4.  A  TWO-DIMENSIONAL  MODEL  WITH  CONSERVED  MASS  AND  MOMENT l 


d=0 . 15 


Time  2 


d=0 . 25 


Time  2 


d=0 . 35 


Time  2 


d=0 . 45 


Time  2 


Figure  4.12:  Sound  wave  density  oscillations,  depicted  in  blue,  measured  from  classical  FHP  lattice-gas  simulations 
carried  out  for  a  range  of  reduced  background  densities,  d  =  0.15  to  d  =  0.65.  Initially,  the  density  field  has  a 
sinusoidal  perturbation  of  size  A  =  256  with  an  amplitude  of  5d  =  0.1  and  with  its  wave  vector  directed  along  the  x- 
axis,  k  =  ( k ,  0).  The  simulation  was  carried  out  on  a  256  x  256  triangular  lattice  for  t  =  3000  time  steps,  and  data  was 
sampled  at  every  other  time  step.  The  amplitude  of  the  sound  wave  decays  exponentially  fast,  5p  =  5p0e'lkCst~k  rt, 
where  the  wave  number  is  k  =  ^  and  the  sound  speed  is  cs  =  —L=.  The  damping  constant,  T,  is  determined  by  a 

A  Tv  2 

numerical  fit  to  the  envelope  of  the  oscillation,  and  the  resulting  exponential  curve  is  overplotted  in  red. 
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The  time  evolution  of  a  sinusoidal  mass  density  perturbation  of  wavelength  A  =  256  in  a  256 
size  system  (with  periodic  boundary  conditions  imposed)  was  simulated  with  a  microscopic  FHP 
lattice  gas.  To  measure  the  damping  constant,  the  simulation  was  run  for  3000  time  steps  and  the 
amplitude  of  the  resulting  compressional  standing  wave  was  record.  The  time  series  data  is  plotted 
in  Figure  4.12  for  a  range  of  background  densities.  The  envelope  of  the  oscillation  is  a  decaying 
exponential  curve.  This  decay  constant  is  —  k2T(p).  In  this  way,  numerical  fitting  allows  us  to 
determine  the  sound  damping  constant,  r(p),  as  a  function  of  mass  density.  The  numerical  results 
measured  from  an  FHP  lattice-gas  simulation  are  compared  with  Boltzmann  predictions  for  a  range 
of  reduced  density  from  d  =  0.15  to  0.75.  The  comparison  is  shown  in  Figure  4.13.  The  agreement 
between  simulation  and  theory  is  good. 
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Figure  4.13:  Comparison  of  theory  versus  simulation  for  the  classical  FHP  lattice  gas.  The  analytical  curve 

(plotted  in  red)  is  calculated  from  a  detailed  analysis  of  the  linearized  Boltzmann  equation.  The  sound  damping 
constant  is  predicted  in  the  long  wavelength  limit  as  is  shown  in  Figure  4.11.  The  numerical  data  (plotted  as  circles) 
is  calculated  from  an  FHP  lattice-gas  simulation  using  the  method  of  relaxation  of  a  mass  density  field  standing  wave 
as  is  shown  in  Figure  4.12.  The  two  methods  of  determining  the  coefficient  of  sound  damping  are  in  good  agreement. 


The  constancy  of  the  sound  speed  was  also  using  the  data  taken  from  the  compressional  wave  test. 
It  appears  that  cs  for  the  microscopic  FHP  lattice  gas  is  independent  of  density.  See  Figure  4.14. 
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Figure  4.14:  Plot  of  the  power  spectrum,  p^Pkt  computed  using  the  discrete  Fourier  transform  of  a  mass  density 
field  standing  wave  for  a  range  of  reduced  background  densities  from  d  =  0.15  to  d  =  0.7.  Measurements  were  taken 
from  classical  FHP  lattice-gas  simulations  on  a  256  x  256  grid.  The  size  of  the  standing  wave  is  A  =  256  and  it  oscillates 
and  decays  over  time  (see  Figure  4.12).  The  peak  in  the  power  spectrum  occurs  at  the  same  location  corresponding 
to  a  density-independent  sound  speed. 


Chapter  5 

Conclusion 


We  have  presented  a  review  of  the  classical  lattice-gas  method  that  included  a  description  of  how 
one  may  analytically  predict  certain  fluid- like  behaviors  in  the  long  wavelength  limit.  We  have 
also  included  a  description  of  two  lattice-gas  models  using  one  and  two-dimensional  lattices  as 
examples.  A  linearized  lattice-Boltzmann  equation  at  the  mesoscopic  scale  is  used  to  calculate  the 
dispersion  relations  for  the  sound  mode  and  shear  and  bulk  viscosity  modes  of  the  dynamical  lattice- 
gas  system.  For  small  wave  numbers,  these  dispersion  relations  are  compared  to  numerical  data 
taken  from  both  large-scale  and  slow  modes  present  in  the  two  lattice-gas  models.  Numerical  and 
analytical  predictions  of  the  sound  damping  constant  and  sound  speed  are  in  agreement  and  diffusive 
ordering  of  damping  times  versus  grid  size  is  observed  in  this  context.  Furthermore,  numerical  and 
analytical  predictions  for  the  shear  and  bulk  viscosities  for  the  two-dimensional  lattice-gas  system 
were  compared  and  are  in  good  agreement  as  well. 

The  numerical  efficiency  and  convergence  properties  of  the  lattice-gas  algorithm  are  not  reported 
in  the  chapters  covering  the  numerical  simulations.  Therefore,  it  is  important  to  state  here  that 
the  computational  efficiency  and  the  order  of  convergence  of  the  lattice-gas  algorithm  implemented 
on  classical  general-purpose  computers  is  much  lower  than  that  achievable  by  computational  fluid 
dynamics  codes  implemented  on  the  same  general-purpose  computers  using  high-level  languages 
based  on  floating-point  representations  of  real  valued  quantities.  A  significant  speed-up  of  the 
lattice-gas  algorithm  can  be  achieved  using  simple  special-purpose  computers  [104,  105].  However, 
the  gains  are  not  sufficient  to  warrant  the  continued  construction  of  these  special-purpose  classical 
computers,  except  perhaps  for  use  in  the  narrowly  defined  application  area  of  the  computational 
fluid  dynamics  of  multispecies  and  multiphase  lattice-gas  systems  and  for  use  in  understanding  the 
extent  to  which  reversible  algorithms  (a  special  case  of  unitary  algorithms)  applied  over  massively 
large  data  sets  can  be  used  for  physical  modeling  purposes. 
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CHAPTER  5.  CONCLUSION 


Appendix  A 

Small  Mach  Number  Expansion  of 
the  Occupancy  Probability 


The  single-particle  distribution  function  has  the  form 

(A.l) 

where  the  natural  log  of  the  fugacity 

In  za  =  ap  +  (3ea  ■  p  +  'yE 

(A.2) 

is  a  linear  combination  of  the  conserved  scalar  quantities,  the  mass  p,the  momentum  component 
ea  ■ p  along  the  lattice  direction  ea,  and  the  energy  if  at  a  lattice  site.  The  real  numbered  coefficients 
a,  p,  and  7  are  free  parameters  that  we  will  determine.  It  is  convenient  to  define  the  momentum 
and  energy  independent  part  of  the  fugacity  as 

20  =  eap. 

(A.3) 

Since  fa{z0)  =  d  is  the  reduced  density,  d  =  -fe,  we  must  set 

1-d 

a 

(A.4) 

This  fixes  the  coefficient  a.  To  fix  the  coefficients  /3  and  7,  we  can  specify  two  moments  of  the 
single-particle  distribution  function  as  constraint  conditions.  We  begin  by  Taylor  expanding  the 
single-particle  distribution  function  f(za)  about  z0 

f(za )  =d  +  f(z0)Sz  +  ^ f"(z0)(Sz 2)  +  • 

(A.5) 

The  derivative  of  /  evaluated  at  zQ  are 

/'(*)=  (z+l)2  ^ 

(A.6) 

and 

/  (z)  —  ^  ^  /  (*>)  —  2 d  , 

so 

f(za)^d[l-dSz  +  d2(Sz)2]  . 

(A.7) 

(A.8) 
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To  determine  Sz,  we  begin  by  writing  the  fugacity  in  series  form 


za  =  z0 


E 


(Pea-p)k 


lk=0 


k\ 


E 


(7  E)k 


J  Lk=0 


k\ 


(A.9) 


In  the  subsonic  limit,  p  <C  me,  keeping  terms  only  to  second  order  in  the  velocity,  the  fugacity 
becomes 

1  +  f3ea  -p  +  \(Pea -p)2  (1  +  ~/E)  +  0(v3). 


za  =  z  o 


since  p  ~  v  and  E  ~  v2 .  Then  to  second  order  in  the  velocity,  the  change  in  za  is 


Sza  =  za  -  z0  = 
and  the  square  of  the  change  is 

( Sza )2  = 


1  —  d 
d 


1. 


Pea-P  +  -(Pea  ■  p)2  +  jE 


0(v3) 


1  -d 


P2  (ea  ■ pf  +  0(v3). 


Inserting  the  expressions  for  Sz  and  (Sz)2  into  the  Taylor  expansion  of  f(za)  we  have 


f(za)  =  d{l-(l-d) 

=  d 


Pea  ■  p  +  -  (Pea  -p)2  +  jE 


+  (1  -  d)2  (ea  ■  p)‘ 


1  -  (1  -  d)  (pia-p  +  'yE)  +  -(1  -  d)(l  -  2d)p2  (ea-p) 


(A.10) 


(A.ll) 


(A.12) 


(A.13) 


We  have  the  freedom  to  choose  the  coefficients  P  and  7  to  parameterized  the  distribution  function 
as  we  see  fit  to  satisfy  any  two  constraints.  Consider  a  parameterization  that  fixes  the  value  of  the 
coefficients  P  and  7  by  using  the  following  moments  for  the  mass  density  and  momentum  density 


B 

p  =  m^fa 

(A.14) 

a= 1 

B 

pv  =  rnc'^T  eafa. 

(A.15) 

a=  1 


The  parameterization  may  be  termed  the  non-Galilean  parametrization.  Constraints  Equations  (A.  14) 
and  (A. 15)  are  typically  used  in  the  formulation  of  classical  lattice  gases.  The  single  particle  distri¬ 
bution  function  using  this  non-Galilean  parameterization  was  first  found  in  the  mid  1980’s  by  the 
US  researchers  Wolfram  and  Hasslacher  and  by  the  French  researchers  Frisch,  d’Humieres,  Lalle- 
mand,  Pomeau,  and  Rivet  [1,  3].  Their  derivation  of  Equation  (A. 20)  is  different  then  the  derivation 
presented  in  this  section;  they  used  only  two  free  coefficients  in  the  expression  for  the  fugacity,  one 
for  the  mass  and  the  other  for  the  momentum;  whereas  we  use  three  free  coefficients.  The  reason 
for  using  only  two  free  parameters  is  that  in  the  standard  single-speed  classical  lattice-gas  construc¬ 
tion,  the  energy  is  degenerate  with  the  mass,  so  it  was  deemed  unnecessary  to  keep  a  separate  free 
coefficient  for  the  energy.  However,  it  is  expedient  to  use  a  free  parameter  for  E.  Using  Equa¬ 
tions  (A. 14)  and  (A. 15)  as  constraint  equations  gives  us  a  non-unity  density-dependent  prefactor  in 
the  convective  term  in  the  hydrodynamic  flow  equation. 

Inserting  Equation  (A.13)  into  Equation  (A. 15),  the  odd  term  in  the  distribution  function  ex¬ 
pansion  survives  the  first  moment  sum  over  lattice  directions;  the  odd  term  is  the  one  linear  in  the 
momentum.  This  fixes  the  value  of  P  to  be 
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so  the  distribution  function  becomes 

fa  =  d 


i  _  o  j 

1  +  Dea  •  p  +  —  y— j  (ea  •  p)2  +  (1  -  d)jE 


(A. 17) 


Inserting  Equation  (A.  17)  into  Equation  (A.  14),  all  the  even  terms  that  survive  the  sum  over  lattice 
directions  must  add  to  zero.  This  fixes  the  value  of  7  as  follows 


D  1  -  2d 
2~  1  -d 


p2  —  (1  —  d)^E  =  0 


or 


"/E  =  D 


1  —  2  d  p2 


(1  -  d)2  2  ' 

Therefore,  the  non-Galilean  parameterized  distribution  function  is 

£>(£>  +  2) 


fa=d 


1  +  DeaiPi  + 


-g{d)QaijPiPj 


where  the  density  dependent  prefactor  g{d)  is  defined 

D  1-2  d 
g{d)  =  D  +  2  1-d 

and  the  traceless  second-rank  tensor  Qa  is  defined 

djj 

D  ' 


Qaij  — 


(A.18) 

(A.19) 

(A. 20) 

(A. 21) 

(A. 22) 


Qa  is  an  isotropic  symmetric  tensor.  This  mass-energy  degeneracy  leads  to  an  anomalous  description 
of  the  lattice-gas  fluid’s  behavior.  Let  us  see  why.  The  second  moment  of  Equation  (A. 20)  gives  the 
momentum  flux  density 

B 

rnc2  ^2  eaie-ajfa  =  PSij  +  gpviVj.  (A. 23) 

a=  1 

The  density-dependent  prefactor  g  appears  in  the  nonlinear  convective  term,  so  this  parametrization 
does  indeed  give  rise  to  non-Galilean  fluid  flow.  The  pressure  in  Equation  (A. 23)  has  a  spurious 
quadratic  velocity  dependence 

P  =  pc2(l-g^]. 


(A. 24) 
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